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I.  INTRODUCTION 

The  advent  of  GPS  has  afforded  the  aerospace  controls  engineer  a  powerful  new 
means  of  controlling  air  vehicles.  Current  guidance  schemes  rely  in  some  part  on 
ground  based  radars,  navigational  aides,  beacons,  localizer  beams,  etc.  Guidance 
and  control  of  the  air  vehicle  have  been  developed  seperately  and  then  combined  in 
a  somewhat  adhoc  process.  Precise  tracking  of  inertial  fixed  trajectories  using  an 
onboard  GPS  integrated  GNC  suite  affords  a  quantum  leap  in  autonomous  flight 
capabilities.  The  greatest  impact  of  the  proposed  technology  is  expected  in  the  area 
of  trajectory  tracking  control  for  autonomous  unmanned  vehicles,  and  automatic 
approach  and  landing  of  manned  vehicles  dicussed  next. 

Control  of  the  commercial  air  traffic  throughout  the  country  continues  to  become 
more  and  more  demanding  as  an  increased  number  of  vehicles  vie  for  limited  access 
to  the  major  commercial  aviation  hubs.  Sophisticated  and  expensive  ground  based 
radar  control  facilities  employ  a  large  number  of  personnel  to  individually  instruct 
the  pilots  of  these  aircraft  on  the  trajectories  they  are  to  fly.  Precise  control  of  the 
aircraft  trajectory  such  as  required  by  an  approach  into  an  airfield  requires  constant 
attention  from  a  ground  based  air  traffic  controller.  Furthermore,  ATC  ability  to 
effectively  control  the  aircraft  is  influenced  by  ground  based  radar  coverage  and  navaid 
equipment  available  and  is  often  negatively  influenced  by  atmospheric  conditions. 
Airfields  with  limited  resources  are  often  unable  to  take-off  and  land  aircraft  requiring 
instrument  departures  and  arrivals. 

Flight  patterns  around  major  aviation  hubs  are,  in  general,  inertia!  based  tra- 
jectories. That  is,  an  aircraft  is  required  to  track  a  certain  path  over  the  ground  while 
adhering  to  a  certian  altitude  schedule,  irrespective  of  air  mass  disturbances.  In  some 


cases,  such  as  on  final  approach  to  land  or  when  the  aerodrome  is  situated  among 
significant  terrain  or  cultural  development,  precise  adherence  to  the  desired  trajectory 
is  crucial  for  flight  safety.  In  any  case,  commanding  an  inertial  trajectory  directly 
utilizing  a  GPS  integrated  GNC  system  could  prove  to  be  more  cost  effective  and 
accurate  than  current  methods.  Furthermore,  such  a  guidance  scheme  could  open  up 
many  more  airports  to  significant  commercial  air  traffic  without  requiring  the  capital 
investment  and  maintenance  of  ground  based  radar  facilities. 

Unmanned  air  vehicles  can  be  a  cost  effective  means  of  power  projection.  Addi- 
tionally, in  some  cases  human  physiological  limits  may  prove  to  be  the  limiting  factor 
in  the  performance  of  an  air  vehicle.  The  precision  delivery  of  munitions  becomes 
of  paramount  importance  as  weapons  and  weapon  delivery  platforms  continue  to  in- 
crease in  cost,  thus  limiting  their  numbers.  All  of  these  concerns  can  be  addressed 
by  autonomous  air  vehicles  utilizing  a  GPS  aided  guidance,  navigation  and  control 
suite. 

All  of  these  applications  have  a  common  thread  running  through  them.  While 
the  particulars  of  the  vehicles  may  vary  considerably,  the  intent  is  to  acheive  au- 
tonomous control  of  their  trajectory.  As  a  proof  of  concept,  this  work  presents  a  new 
design  process  for  the  synthesis  of  a  guidance,  navigation,  and  control  system  for  a 
UAV  named  Bluebird.  The  function  of  the  this  GNC  system  is  to  track  inertial  tra- 
jectories. Bluebird  is  a  UAV  operated  at  the  Unmanned  Air  Vehicle  Lab  at  the  Naval 
Postgraduate  School.  It  has  a  12.5  foot  wingspan  and  a  20  pound  payload  capability, 
and  is  currently  being  equipped  with  a  full  avionics  suite,  including  IMU,  GPS,  and 
air  data  sensors. 

The  design  process  began  with  the  development  of  a  nonlinear  dynamic  model  of 
Bluebird  implemented  in  SIMULINK.  A  typical  cruise  flight  condition  was  chosen  as 
the  point  for  linearization.  After  linearization  of  the  nonlinear  model,  the  work  cen- 


tered  on  the  design  of  a  linear  controller.  LQR  (Linear  Quadratic  Regulator)  synthesis 
approach  was  used  since  it  provides  an  intuitive  means  of  synthesizing  a  multivariable 
controller  within  the  framework  of  real  world  design  constraints.  Following  the  de- 
sign of  the  LQR  controller,  the  challenge  of  implementing  the  linear  controller  on  the 
nonlinear  plant  was  addressed.  A  novel  method  of  converting  commands  and  outputs 
from  inertial  to  body  reference  cooridinates  was  used  [Ref.  10].  This  method  achieves 
automatic  recruiting  of  the  actuators,  while  preserving  a  certain  linearization  prop- 
erty. Next,  the  accuracy  of  the  nonlinear  simulation  was  enhanced  with  the  addition 
of  high  fidelity  models  of  sensors  used  onboard  Bluebird.  Additionally,  Kalman  filters 
were  designed  in  order  to  provide  optimal  state  estimates.  Finally,  the  performance 
of  the  controller  was  evaluated  in  simulations  with  the  full  nonlinear  model. 


II.  DEVELOPMENT  OF  THE  DYNAMIC 

MODEL 

The  development  of  an  integrated  guidance,  navigation,  and  control  system 
required  a  high  fidelity  nonlinear  model  of  the  aircraft  dynamics.  The  discussion 
begins  with  an  explanation  of  nomenclature,  abbreviations,  and  a  definition  of  frames 
of  reference. 

A.     REFERENCE  FRAMES 

Three  different  reference  frames  are  used  in  this  report.  They  are: 

•  Local  Tangent  Plane  or  Inertial  Reference  Frame 

•  Body- Fixed  Reference  Frame 

•  Wind  or  Flight  Path  Reference  Frame 

1.     Local  Tangent  Plane  Reference  Frame 

The  position  of  the  air  vehicle  must  be  maintained  with  respect  to  the  local 
tangent  plane  coordinate  system.  This  coordinate  system  is  formed  by  extending  a 
ray  from  the  center  of  the  earth  to  its  surface.  A  plane  is  attached  tangent  to  the  point 
of  intersection  of  the  ray  with  the  Earth's  surface.  While  it  is  somewhat  arbitrary,  for 
our  purposes  here  it  will  be  convenient  to  define  the  positive  x  direction  as  pointing 
east,  the  positive  y  direction  as  pointing  north,  and  the  positive  z  direction  as  pointing 
up.  This  is  depicted  in  Figure  2.1. 

For  the  purposes  of  this  development,  the  rotation  of  the  earth  and  its 
associated  Coriolis'  forces  can  be  ignored  and  the  local  tangent  plane  reference  frame 
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Figure  2.1:  Local  Tangent  Plane  Coordinate  System 

can  be  considered  to  be  an  inertia!  reference  frame.    In  this  work,  {/}  is  used  to 
represent  the  inertia!  reference  frame. 

2.     Body-Fixed  Reference  Frame 

The  body-fixed  reference  frame  is  a  right  hand  orthogonal  system  with 
the  origin  at  the  center  of  gravity  of  the  air  vehicle.  The  positive  x  direction  points 
towards  the  nose.  The  positive  y  direction  points  out  the  right  wing  and  the  positive  z 
direction  points  towards  the  bottom  of  the  air  vehicle.  The  velocity  of  the  air  vehicle 
with  respect  to  the  inertia!  reference  frame,  resolved  along  the  x,  y,  and  z  axis  of  the 
body-fixed  reference  frame,  are  termed  u,  v,  and  w,  respectively.  The  angular  rate  of 
rotation  of  the  air  vehicle  with  respect  to  the  inertia!  reference  frame,  resolved  in  the 
body-fixed  reference  frame,  are  called  p,  q,  and  r,  respectively.  Positive  values  for  the 
forces,  moments,  angular  rates,  and  linear  velocities  in  the  body-fixed  reference  frame 
are  shown  in  Figure  2.2.  The  abrieviation,  {#},  is  used  to  represent  the  body-fixed 
reference  frame. 
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Figure  2.2:  Body-Fixed  Coordinate  System  [Ref.  11] 

3.     Flight  Path  or  Wind  Reference  Frame 

The  wind  reference  frame  is  also  a  right  hand  orthogonal  system  with  its 
origin  at  the  center  of  gravity,  e.g.,  of  the  air  vehicle.  The  x  axis  is  aligned  with  the 
velocity  vector  of  the  air  vehicle.  The  orientation  of  the  wind  reference  frame  with 
respect  to  the  body-fixed  reference  frame  is  defined  in  terms  of  the  angles  a  and  /?. 
The  equations  for  a  and  (3  are  given  below. 


a  =  tan     (w/u) 


(2.1) 


and 


/3  =  sin-l(v/V) 


(2.2) 


where  the  vectors  u,  v,  w,  and  V  are  velocity  components  of  the  air  vehicle  defined  in 
Figure  2.3.  The  abrieviation,  {W},  is  used  to  represent  the  wind  reference  frame. 


Figure  2.3:  Wind  or  Flight  Path  Reference  Frame  [Ref.  11] 


B.     COORDINATE  TRANSFORMATIONS 

In  order  to  use  these  three  coordinate  systems,  one  must  be  able  to  transform 
between  them  freely.  The  Euler  angles,  $,0,  and  #,  termed  roll,  pitch,  and  yaw, 
are  defined  in  order  to  express  the  orientation  of  the  body-fixed  reference  frame  with 
respect  to  the  inertial  reference  frame.  For  the  purposes  of  this  development,  a  3-2-1 
Euler  angle  transformation  will  suffice  as  its  singularity  occurs  at  0  equal  to  90  de- 
grees. The  3-2-1  transformation  is  given  without  explanation  but  a  good  development 
of  Euler  angle  transformations  in  general  can  be  found  in  [Ref.  6].  The  nature  of 
the  angular  rotation  is  more  apparent  when  the  transformation  is  expressed  as  the 


product  of  three  rotation  matricies.  In  the  case  of  a  3-2-1  rotation  sequence,  the  three 
matrices  in  Equation  2.3  correspond  to  rotations  about  the  yaw,  pitch,  and  roll  axes 
of  the  air  vehicle.  Of  course,  the  three  matrices  can  be  multiplied  out  for  an  analytic 
result  contained  in  a  single  matrix  ,  although  the  resulting  matrix  is  somewhat  busy 
to  inspect.  In  any  case,  the  transformation  between  a  free  vector  resolved  in  the  in- 
ertial  reference  frame  and  the  same  vector  resolved  in  the  body-fixed  reference  frame 
is  given  by: 


lV  = 


cos  \P      sin  $    0 

-  sin  $    cos  #    0 

0  0        1 


cos  0    0    —  sin  0 

0       1         0 
sin  0    0      cos  0 


1         0  0 

0      cos  $      sin  $ 
0    —  sin  $    cos  $ 


V    (2.3) 


where  lV  is  a  free  vector  resolved  in  {/}  and  BV  is  the  same  vector  resolved  in  {B}. 
The  inverse  is  also  denned.  Conveniently,  since  the  transformation  is  orthonormal, 
the  inverse  is  simply  the  transpose  of  the  rotation  matrices  shown  in  Equation  2.3. 

Not  all,  transformations,  however,  are  orthonormal.  Of  particular  interest  is 
the  case  of  angular  rotation  rates.  The  body-fixed  reference  frame's  angular  rate  of 
rotation  with  respect  to  the  inertia!  reference  frame  can  be  related  to  the  rate  of 
change  of  the  Euler  angles  by  a  transformation  matrix.  The  development  is  straight 
forward  and  is  fully  explained  in  [Ref.  11].  The  final  transformation  matrix  from  p, 
q,  r  to  the  time  rate  of  change  of  the  Euler  angles,  $, 0,  ^,  is  given  by: 


0 


1    sin  $  tan  0    cos  $  tan  0 
0         cos  $  —  sin  $ 

0    sin$  sec0     cos$  sec0 


P 

q 

r 


(2.4) 


By  integrating  Equation  2.4,  the  time  history  of  the  Euler  angles  can  be  obtained. 

Aerodynamic  forces  and  moments  are  often  calculated  using  stability  and  control 
derivatives  defined  with  respect  to  the  wind  reference  frame.  The  angles,  a  and  /?, 
define  the  orientation  of  the  wind  reference  frame  to  the  body-fixed  reference  frame. 
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Therefore,  a  transformation  matrix  can  be  obtained  that  relates  a  free  vector,  such  as 
lift  or  drag,  resolved  in  {W}  to  the  same  vector  resolved  in  {B}.  The  transformation 
is  expressed  as: 


v  = 


cos  a  cos  ft    —  cos  a  sin  /?    —  sin  a 

sin  f3  cos  /?  0 

sin  a  cos  (3    —  sin  a  sin  /?      cos  a 


w 


V  (2.5) 


where  BV  is  a  free  vector  resolved  in  {£}  and  ^V  is  the  same  vector  resolved  in 
{W}. 

C.     NOTATION 

Some  standardized  abbreviations  will  simplify  the  development  of  the  nonlinear 
kinematic  model  of  the  air  vehicle.  This  short-hand  is  used  in  the  field  of  robotics 
where  multiple  frames  of  reference  are  common  [Ref.  5]. 

•  Peg  represents  the  position  vector  from  the  origin  of  the  local  tangent  plane  to 
the  center  of  gravity  of  the  air  vehicle. 

•  Vcg  and  Bacg  represent  the  velocity  and  acceleration,  measured  at  the  center  of 
gravity  of  the  air  vehicle,  with  respect  to  {/},  resolved  in  {B}.  The  components 
of    Vcg  are  commonly  termed  u,  v,  and  w. 

•  Ivcg  and  Iacg  represent  the  velocity  and  acceleration,  measured  at  the  center  of 
gravity  of  the  air  vehicle,  with  respect  to  {/},  resolved  in  {/}. 

•  bub  is  the  angular  velocity  of  the  {B}  coordinate  system  with  respect  to  {/}, 
resolved  in  {B}.  The  components  of  bub  are  commonly  termed  p,  q,  and  r. 

•  iub  represents  the  angular  velocity  of  the  {B}  coordinate  system  with  respect 
to  {/},  resolved  in  {/}. 


•  gR  represents  the  transformation  matrix  used  to  express  a  free  vector  resolved 
in  {B},  in  {/}.  The  inverse  is  represented  by  BR 

•  yyR  represents  the  transformation  matrix  used  to  express  a  free  vector  resolved 
in  {W},  in  {B}.  The  inverse  is  represented  by  g^R. 

•  BF  and  BN  denote  the  total  external  inertial  force  and  moment  acting  on  the 
body  resolved  in  {B}. 

•  *F  and  *N  denote  the  total  external  inertial  force  and  moment  acting  on  the 
body  resolved  in  {/}. 

•  BL  is  the  inertial  angular  momentum  of  the  body  resolved  in  {B}. 

•  IL  is  the  inertial  angular  momentum  of  the  body  resolved  in  { /}. 

•  Given  a  vector  u,  its  derivative  with  respect  to  {B}  is  denoted  as  Jj(u) 
and  its  derivative  with  respect  to  {/}  is  denoted  as  (v) 

D.     RIGID  BODY  EQUATIONS  OF  MOTION 

In  general,  an  avionics  suite  on  a  modern  air  vehicle  utilizes  a  strapdown  IMU. 
A  strapdown  IMU,  as  the  name  implies,  maintains  a  constant  orientation  in  the  body- 
fixed  reference  frame.  The  output  of  the  sensors  on  the  IMU  are  resolved  in  {B}. 
Therefore,  among  other  reasons,  it  is  most  convenient  to  develop  the  equations  of 
motion  of  the  air  vehicle  in  the  body-fixed  reference  frame. 

1.     Linear  Motion 

An  application  of  Newton's  Law  to  linear  motion  of  a  body  states  that 
the  total  external  force  applied  to  a  body  is  equal  to  the  mass  of  the  body  times  its 
inertial  acceleration.  This  could  be  written  in  the  inertial  reference  frame  as: 
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where 

'a    =    '*»,  (2.6) 

or  in  the  body-fixed  reference  frame  as  follows: 


B  j?     _     __  B„ 
r     =    m    a 

=    mfRti) 


eg 


=    m  **<.,.  (2.7) 

Coriolis'  theorem  can  be  used  to  relate  the  inert ial  and  body  accelerations  of  the  air 
vehicle  as  follows: 

\  =  -\  +  BWflx\S)  (2.8) 

where  the  difference  in  the  derivatives  is  explained  in  Chapter  II,  part  C.  Equation  2.8 
can  be  substitued  into  Equation  2.7  in  order  to  obtain  the  desired  expression  for  the 
sum  of  the  external  inertia!  forces  resolved  in  the  body-fixed  reference  frame. 


BF    =    m{-Bvcg+BuBxBvcg) 

=    m  -JlBvc9  +  rn  (BU)B  x  Bvcg).  (2.9) 

at 

2.     Angular  Rotation 

Euler's  law  for  the  conservation  of  angular  momentum  at  the  center  of 
gravity  states  that: 
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'Leg  =  'Nc,,,  (2.10) 

where  ILcg  is  the  angular  momentum  of  the  air  vehicle  with  respect  to  {/}  and  JNcg 
is  the  total  moment  applied  to  the  air  vehicle.  Equation  2.10  can  be  written  in  the 
body-fixed  reference  frame  as: 

BLcg  =  fR'Nc,.  (2.11) 

Coriolis'  theorem  can  be  used  again  to  expand  BLcg  obtaining: 

BLeg  =  ^BLcg+BuBxBLcg.  (2.12) 

It  can  be  shown  that  the  angular  momentum,  BLcg,  of  the  air  vehicle  is  the  product 
of  an  inertia  tensor,  denned  as  JB,  and  the  body's  angular  velocity,  BuB,  where  we 
ignored  all  spining  elements.  Substituting  this  definition  of  BLcg  into  Equation  2.12, 
results  in: 

BLcg  =  jt(JBBUB)  +  %J  x  Jbbujb.  (2.13) 

Recall  that  B L^  =  f&Ncg  =B  N^.  Using  this  relationship,  Equation  2.13  can  be 
equivalently  expressed  as: 

BNcg  =  ^-(JBBuB)  +  BuB  x  JBBu>B  (2.14) 

at 

3.     External  Forces  and  Moments 

Equation  2.9  and  Equation  2.14  from  the  preceding  sections  can  be  com- 
pactly expressed  as  follows. 
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Bp 


m 


ftBvcg  +  m  (BuB  x  Bvcg) 


B 


JBTtUUB  +  UU)B  X  JBBUB 


B, 


dt 


(2.15) 


By  bringing  derivative  terms  in  Equation  2.15  to  the  left  hand  side,  we  obtain, 


d_ 
dt 


v, 


eg 

3ub 


-bujb  X  Bvcg 


+ 


BF 


-Jg1  BuB  x  Jbbujb    +    JE^8 


N 


(2.16) 


Next,  the  forces  and  moments  in  Equation  2.16  can  be  expanded  as  follows: 


Bp 

BN 


FgRAVITATIONAL  +      FpBOPULSIVE  +  BF AERODYNAMIC 

Npropulsive  +  bNa 


'AERODYNAMIC 


,      (2.17) 


where 


FgRAVITATIONAL 
PROPULSIVE 
AERODYNAMIC 
^PROPULSIVE 
AERODYNAMIC 


bFf 

BFJ 


BN 


force  due  to  gravity 

force  due  to  engine's  thrust 

force  generated  by  aerodynamic  surfaces 

moment  generated  by  engine's  thrust 

moment  generated  by  aerodynamic  surfaces 


The  gravitational  force  expressed  in  {/}  is  given  by 

Fgravity  = 
where  "g"  is  the  gravitational  constant.  Then 


0 

0 

mg 


B 


Bnl 


Fgravity  =  j  R  Fgravity- 


(2.18) 


The  expansion  of  the  propulsive  forces  and  moments  is  simplified  by  consid- 
ering the  case  of  centerline  thrust.  In  that  case,  no  external  moments  are  generated, 
i.e.,     Npropulsive  =  0,  and  the  propulsive  forces  can  be  expressed  as: 


B 


Fprop  = 


T 
0 
0 


(2.19) 
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where  T  represents  the  thrust  of  the  powerplant. 

Aerodynamic  forces  and  moments  are  commonly  calculated  using  nondi- 
mensional  stability  and  control  derivatives.  These  derivatives  are  obtained  by  approx- 
imating the  aerodynamic  forces  and  moments  acting  on  the  air  vehicle  using  a  Taylor 
series  expansion  about  a  given  trim  point.  Typically,  values  for  these  derivatives  are 
available  for  the  first  order  terms  of  the  expansion  only.  Sometimes,  a  few  second 
order  terms  are  available,  such  as  the  terms  associated  with  a  and  /?  [Ref.  7].  All 
other  higher  order  terms  are  usually  ignored  in  this  approximation. 

In  general,  the  aerodynamic  forces  and  moments  acting  on  the  air  vehicle 
are  computed  as  follows.  The  nondimensional  stability  and  control  derivatives  are 
dimensionalized  by  multiplication  by  the  appropriate  constants,  such  as  wing  span, 
dynamic  pressure,  chord,  etc.  The  dimensional  derivative  is  then  multiplied  by  the 
perturbations  of  each  aerodynamic  variable  or  control  deflection  from  its  nominal 
trim  point.  The  summation  of  the  forces  and  moments  due  to  all  of  the  aerodynamic 
variables  and  control  deflections,  in  addition  to  the  trim  value  of  the  forces  and 
moments,  results  in  the  total  aerodynamic  force  and  moment  acting  on  the  air  vehicle. 
4.     The  State  Space  Representation 

In  order  to  implement  the  dynamic  model  in  a  state-space  form  suitable  for 
numerical  simulation,  the  states  of  the  model  need  to  be  chosen.  This  is  somewhat 
arbitrary  and  many  choices  will  work  so  long  as  consistency  is  maintained  in  the 
approach. 

As  was  evident  from  the  development  of  the  rigid  body  equations  of  motion, 
the  body-fixed  reference  frame  is  the  most  convenient  coordinate  system  in  which  to 
define  the  states.  The  first  three  states  are  defined  as  the  inertial  velocity  of  the  air 
vehicle  resolved  in  the  body-fixed  reference  frame.  These  are  abrieviated  as  u,  v,  and 
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w  or  more  compactly  as  Bvcg.    The  fourth  through  sixth  states  are  defined  as  the 
angular  velocity  of  the  air  vehicle  with  respect  to  {/}  resolved  in  {B}.    These  are 
abrieviated  as  p,  q,  and  r  or  more  compactly  as  bub  ■ 
Control  inputs  are  represented  by  the  vector  A: 


A  =  [6e,  6r,  6a] 


(2.20) 


where  6C,  6r,  and  6a  are  the  elevator,  rudder,  and  aileron  inputs,  respectively. 
Control  inputs  for  the  throttle  are  represented  by  6trt> 

Typically,  the  terms  in  the  Taylor  series  expansion  for  the  aerodynamic 
stability  derivatives  are  partial  derivatives  with  respect  to  u/U,  a,  /?,  p,  q,  and  r, 
where  U  is  the  magnitude  of  the  air  vehicle's  velocity  vector  and  a  and  /3  are  the 
angles  defining  the  orientation  of  {B}  with  respect  to  {W}  [Ref.  3].  The  last  three 
variables,  p,  qr,  and  r,  are  states  of  the  model.  The  first  three  can  be  represented  as 
a  combination  of  states  in  the  model  as  follows.  First  note  that, 


U 
0 
0 


Wr>B 


_     w  no,, 
—    B  "•    vcgi 


(2.21) 


and  for  small  values  of  angles  a  and  /?,  a  is  approximately  equal  to  w/U  and  /?  is 
approximately  equal  to  v/U. 

It  turns  out,  the  stability  derivatives  can  be  placed  in  matrix  form  as 
follows: 


dC 
dx' 


Clv  Clp  Cl0  Clp  CLq  Clt 

Cyu  Cyp  Cya  Cyp  CyCi  Cyt 

Cdu  Cdp  Coa  Cdp  Cd„  Cdf 

Ciu  Ch  Ci«  ^ip  Cf«  C'r 

^mu  ^rnp  ^ma  ^mp  ^m,  ^mp 

Cnu  Cnp  Cna  Gnp  Cn,  Onr 


Similarily,  for  the  control  derivatives: 
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\Cl„ 

Cl,t 

CL,a 

Cy,. 

Cvlr 

Cy,. 

dC 

Co,. 

CDir 

Cota 

dA~ 

c«. 

Cilr 

CiSa 

Cmge 

Crnsr 

CrnSa 

.  cnge 

Cntr 

Cnfa    . 

where  x'  is  the  vector  composed  of  u/U,  a,  /?,  p,  q,  and    r.  Now  the  aerodynamic 
forces  and  moments  can  be  expressed  as  follows: 


BN 


Faebo 

AERO 


=  qS 


%R     0 
0      frR 


^+SM'*+l>i+lfA>'  <2-22> 


where  M\  q,  and  5  are  matrices  used  to  dimensionalize  the  stability  and  control 
derivatives  and  convert  the  state  vector,  x,  to  x': 


=  I  u    v 


w    p    q    r 


]T. 


and 


q  =  dynamic  pressure, 
5  =  diag{— 5,  5,  —5, 56,  Sc,  56}, 


x'  =  M'x, 


M'  =  diag{l/FT,  1/Vt,  1/Vt,  b/2VT,  c/2VT,  b/2VT}, 


x'  =  M'x, 

M'  =  diag{0,C/(2VT),6/(2yT),  0,0,0}. 

Equation  2.16  can  now  be  further  expanded  using  expressions  of  the  forces 
and  moments  derived  above. 
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d_ 
dt 


B 


'eg 
UB 


B 


0 


0 
-bJb\b»b*bJbbub) 

0    BJsl 


Vcg 

BF 

BN 


+ 


(2.23) 


where 


{ 


BF 

BN 


-{ 


BF 


GRAV 
0 


+ 


r  bf 


PROP 
0 


Strt+ 


&R     0 

0  ^R 


qS  {  CFO  +  %M'x  +  %M>x  +  |f  A  }  }  }  .  (2.24) 


Notice  that  there  is  a  state  derivative  term  on  the  right  hand  side  of  Equa- 
tion 2.23  due  to  the  second  order  terms  in  the  Taylor  series  expansion  of  the  aerody- 
namic forces  and  moments  in  Equation  2.24.  By  bringing  it  to  the  left  hand  side  of 
Equation  2.24  and  combining  terms,  we  obtain: 


d_ 
dt 


Ucg 

*UB 


-'"{[I 


-bubx 


MJ^TqS^M' 

bFprop 
0 


Bx 

B 


'eg 
U>B 


0 
-*Jb1(bubXbJbbub) 

bFgrav 


+  M-{ 


0 


+ 
+ 


Stn  +  ^TqS{CFo  +  &A) 


(2.25) 


where: 


B  rp  _ 

w1  — 


ByR      0 
0      &R 


and  Mi  = 


and 


X  =  /e-M,-Wgf 


m       0 
0     BJB 


AT. 


(2.26) 


Equation  2.25  expresses  the  derivative  of  the  first  six  states  of  the  nonlinear  model 
in  matrix  form.  It  is  solved  in  the  user  denned  MATLAB  Fen  block,  state-deriv.m, 
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which  is  available  for  inspection  in  Appendix  A.  The  physical  parameters  specific 
to  Bluebird  are  stored  in  a  MATLAB  .m  file,  Bluebird-data. m,  and  are  called  from 
within  the  function,  state-deriv.m  .  In  order  to  change  the  vehicle  being  modeled, 
one  simply  changes  the  constants  in  Bluebird-data.m  . 

Next,  the  Euler  angles  were  added  as  the  additional  three  states  of  the 
model.  The  time  history  of  the  Euler  angles  is  denned  in  Equation  2.4  and  written 
in  compact  form  as: 


A  =  S(A)%,, 


(2.27) 


where 


A  = 


0 


and 


5(A)  = 


1    sin  $  tan  0    cos  $  tan  0 
0         cos  $  —  sin  $ 

0    sin$  sec0     cos$  sec0 


The  user  defined  MATLAB  Fen  block,  eu/.m,  solves  Equation  2.27,  the  details  of 
which  can  be  found  in  Appendix  A. 

Finally,  the  inertia!  position  of  Bluebird  can  be  computed  as  follows. 


lb    _/ „      -ipB„ 

rcg  —     vcg  —  B*l    V 


eg- 


(2.28) 


The  rotation  matrix,  JBR,  is  implemented  in  a  MATLAB  Fen  block, 
posb2i.m,  in  order  to  convert  the  vector  Bvcg  to  the  vector  Ivcg.  The  vector  JVcg 
is  then  integrated  to  obtain  a  time  history  of  the  position  of  the  air  vehicle  in  the 
local  tangent  plane.  To  increase  fidelity  of  the  simulation,  a  first  order  model  of  the 
four  actuators,  levator,  Kidder,  ^aileron,  ^throttle,  is  included  with  a  time  constant 
of  1/12.    The  resulting  nonlinear  dynamic  model  now  has  sixteen  states  which  are 
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computed  using  Equations  2.25,  2.27,and  2.28;  summarized  here  for  clarity: 


states  =[u    v    w    p    q    r    $    Q    V    Xpos    Ypos    Zpos    6epo,    6Tpot    6apot    6tpo, 

=  [X    B"l    AT    PT    Ad*. 


il 


poa 


d_ 
dt 


B 


'eg 
U>B 


»-{[[ 


-flo;sx 


0 


-bJb\bub*bJbbub) 


Mj^TqS^M'  ] 


B 


'eg 
UB 


+  Mr1  { 


Fgrav 
0 


+ 


+ 


BFP 


ROP 


0 


Strt  +  ^TqS{CFo  +  ff  A) 


}}■ 


(2.29) 


A  =  S(A)BuB, 


(2.30) 


eg- 


(2.31) 


The  SIMULINK  diagram  of  the  nonlinear  model  is  shown  in  Figure  2.4. 
5.     Trim  and  Linearization 

For  linear  controller  design,  the  nonlinear  equations  above  must  be  trimmed 
and  linearized  for  a  typical  cruise  flight  condition  for  Bluebird.  A  SIMULINK  tool  is 
available  which  can  be  used  to  find  equilibrium  points  of  nonlinear  dynamic  models. 
The  user  specifies  which  states  and  control  inputs  are  to  remain  fixed  along  with 
their  stationary  values  and  the  trim  routine  searches  for  values  of  the  state  and  input 
vector  for  which  the  derivative  of  the  state  vector  equals  zero.  With  the  trim  condition 
known,  another  SIMULINK  tool  perturbs  the  states  around  the  specified  trim  point 
in  order  to  find  the  rate  of  change  of  the  states  and  control  inputs  (Jacobian).  The 
resulting  linear  model  is  returned  in  state  space  format.  Since  Equation  2.28  can  only 
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t- 
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{B}to{l} 

First  Order 
Actuators 

Figure  2.4:  SIMULINK  Nonlinear  Sixteen  State  Dynamic  Model  of  an  Air 
Vehicle 

be  trimmed  for  v^  =  0,  not  a  typical  flight  condition,  we  will  trim  Equations  2.25 
and  2.27,  and  then  include  Equation  2.28  for  linearization. 

We  are  interested  in  trimming  the  model  in  velocity.  Hopefully,  the  deriva- 
tive of  the  position  states  will  never  equal  zero  in  flight.  Also,  in  trim,  the  control 
inputs  are  equal  to  the  actuator  positions.  Therefore,  actuator  states  are  also  removed 
from  the  nonlinear  model.  The  nine  state  nonlinear  model  of  Bluebird  used  for  trim 
is  shown  in  Figure  2.5. 

A  typical  cruise  flight  condition  for  Bluebird  is  given  by: 
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Mux 


d(Euler)/dt=S(Euler)*w 
Equation  1.25 


u 

V 

w 

p 
q 

r 
phi 

theta 
psi 


Figure  2.5:   SIMULINK  Nonlinear  Nine  State  Dynamic  Model  of  an  Air 
Vehicle 

•  flight  speed  equal  to  73  feet  per  second 

•  flight  path  angle  equal  to  zero 

•  wings  level  attitude 

The  nonlinear  model  depicted  in  Figure  2.4  was  trimmed  at  this  condition. 
The  trim  values  of  the  nine  states,  u,  v,  w,  p,  q,  r,  $,  0,  and  #,  and  the  four 
control  inputs,  £e,  6r,  6a,  and  6t  were  returned.  While  the  sixteen  state  nonlinear 
model  cannot  be  trimmed  in  position,  it  can  be  linearized  at  an  arbitrary  position. 
The  origin  of  {/}  is  conveniently  chosen.  These  values  were  then  used  to  linearize 
the  complete  nonlinear  model,  including  position  and  actuator  states.  The  resulting 
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linear  model  of  Bluebird  and  numerical  values  for  the  trim  condition  are  included  in 
the  Appendix  B. 
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III.  THE  LINEAR  QUADRATIC  REGULATOR 

DESIGN 

The  previous  chapter  developed  a  nonlinear  model  of  Bluebird,  which  was  used 
to  derive  a  linear  model  for  a  cruise  flight  condition.  In  this  chapter,  a  linear  dynamic 
controller  is  developed  to  provide  trajectory  tracking  for  the  linear  model.  LQR 
methodology  was  selected  to  design  the  controller.  Based  on  design  requirements,  an 
intuitive  means  of  manipulating  the  LQR  gains  is  presented.  See  [Ref.  9]  for  details. 
The  following  is  a  brief  review  of  the  properties  of  an  LQR  controller  utilized  in  this 
design  process. 

A.     LQR  OVERVIEW 

Consider  the  linear  system 

x    =    Ax  +  Bu 
g  =\  z    =    CiX  +  DlU  (3.1) 

i  y    =    x 

where  x  €  #n,  u  €  BT  and  z  £  B?. 

Suppose  C\D\  =  0  and  D\  is  full  column  rank.  Then, 

zTz  =  xTC?ClX  +  uTD[Dxu. 
Define  a  cost,  J,  as  follows: 

J  =  J  (zTz)dt  =  J  {xTC?ClX  +  uTDlDxu)dt  (3.2) 

and  let  Q  =  CfCi,  and  R  =  D\DX.  Note:  Q  >0  and  R  >  0. 

Assume  (C\,A)  is  observable  and  (A,  B)  is  controllable.  Consider  Figure  3.1. 
The  standard  LQR  problem  is  to  find  a  controller,  u  =  K(s)x,  such  that  the  feed 
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back  system  in  Figure  3.1  is  internally  stable  and  J  is  minimized.  It  turns  out,  one 
such  controller  uses  a  constant  gain,  u  =  —Kx,  where 


K  =  -R-lBJP, 
and  P  solves  the  Alegbraic  Riccatti  Equation: 


-iDr 


A1  P  +  PA  -  PBR^B1  P  +  Q  =  0 


(3.3) 


(3.4) 


Figure  3.1  shows  the  feedback  interconnection  of  the  plant  Q  and  the  controller  K. 
Here  the  inputs  w  can  include  commands  and  disturbances. 

Q  ci 


w 

Q 

t  z 

u 

X 

K 

Figure  3.1:  Standard  LQR  feedback  configuration. 

It  turns  out  that  the  controller,  K,  has  guaranteed  simultaneous  phase  and  gain 
margins  of  no  less  than  60  degrees  and  3dB,  respectively  [Ref.  12].  Furthermore,  the 
controller  has  asymptotic  properties  which  are  exploited  in  the  design  process  and 
are  discussed  next. 
Define  the  Hamiltonian  matrix,  H ,  as: 


H  = 


Let  T  be  given  by: 


\a  -br-xbt 

Q        -AT 
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T  = 


then: 


r-i  = 


/     0 

p  I 

I      0 

-p  I 

where  P  solves  the  Riccatti  equation,  Equation  3.4. 
Note, 


A-BR~lBTP        -BET1BTP 

0  -AT  +  PBR~lBT 


T-iHT  _ 
Therefore  the  eigenvalues  of  H  are  the  roots  of  the  following  polynomials: 


det{sI-H)  =  det{sI-T-xHT)  =  det(sI-A-\-BR-1BTP)det((sI+AT-PBR-1BT) 

Clearly,  the  eigenvalues  of  the  Hamiltonian  consist  of  the  eigenvalues  of  Q  d  and  their 
unstable  reflections  about  the  imaginary  axis.  Let  R  =  pR\,  R\  >  0,  then; 


si  -A    (l/p)R^1BT 
0  sl  +  A 


det(sl  -H)  = 
It  can  be  shown  that  the  det(sl  —  H)  can  be  equivalently  expresses  as  [Ref.  9]: 


det{sI-H)  =  -lndet{sI-A)det{-sI-A)det{I+{l/p)R;1BT{-sI-AT)-1C{C1(sI-A)-1B) 


Let 


and 


Then 


0(a)  =  det(sl  -  A) 


0{s)  =  Cx(aI-A)-lB. 


det(sl  -H)  =  -r<f>(s)<f>(-s)det(I  +  {l/p)R^ l6{-s)6{s))) 


(3.5) 
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The  SISO  example  will  best  demonstrate  what  occurs  to  the  eigenvalues  of  Q  d 
as  p  is  varied  from  0  to  oo.  Let 

e(s)  =  *M/*w 

where  tp(s)  are  the  zeros  of  0(s).  Then, 

det(sl  -H)  =  -l»tf  *)*(-«)(/  +  (l/pW*)1>(-')/<K*)<K-*))' 
It  follows  that  the  eigenvalues  of  H  are  the  roots  of, 

<f>(s)<f>(-s)  +  (l/pMsM-s).  (3.6) 

Consider  the  feedback  system  shown  in  Figure  3.2.  This  is  standard  configuration  for 
SISO  root  locus  analysis  and  its  characteristic  equation  is: 

<p(s)<f>(s)  +  (l/p)1>(s)1>(-3), 

exactly  the  same  as  Equation  3.6.   Since  we  know  that  the  stable  eigenvalues  of  H 
are  the  eigenvalues  of  Q  c/,  standard  root  locus  techniques  show  that 
as  p  goes  to  0  [Ref.  9]: 

•  p  eigenvalues  of  Q  <j  go  to  the  stable  zeros  of  C\{sl  —  A)~XB  or  the  stable 
reflection  of  the  unstable  zeros  of  C\{sl  —  A)~XB,  where  p  is  the  number  of 
zeros  of  Cx{sl  -  A)~lB. 

•  n-p  eigenvalues  of  Q  d  go  to  -oo  in  Butterworth  patterns. 
as  p  goes  to  oo 

•  n  eigenvalues  of  Q  d  go  to  the  stable  eigenvalues  of  A  and  the  stable  reflections 
of  the  unstable  eigenvalues  of  A. 
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Figure  3.2:  Feedback  Configuration  for  Root  Locus  Analysis. 

B.     DESIGN  REQUIREMENTS 

These  properties  of  an  LQR  controller  were  utilized  in  a  design  process  in  order 
to  meet  the  following  design  requirements. 

1.  Zero  Steady  State  Error 


•  Achieve  zero  steady  state  values  for  all  error  variables  in  response  to  ramp 
commands  in  position  along  the  x,  y,  and  z  inertial  axes.  Note  that  a  ramp 
command  in  postion  corresponds  to  a  constant  heading,  constant  velocity 
trajectory. 

2.  Bandwidth  Requirements 

•  The  input-output  command  response  bandwidth  (command-loop  band- 
width) along  any  of  the  three  command  channels  should  be  no  greater 
than  1  radian  per  second  and  no  less  than  1/10  radian  per  second. 

•  The  control-loop  bandwidth  should  not  exceed  12  radians  per  second  for 
the  elevator  and  aileron  actuators,  and  5  radians  per  second  for  the  throttle 
actuator.    These  numbers  represent  80  %  of  the  corresponding  actuator 
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bandwidths  and  shall  ensure  the  actuators  are  not  driven  beyond  their 
linear  operating  range. 

3.  Closed  Loop  Damping 

•  The  dominant  closed  loop  eigenvalues  should  have  a  damping  ratio  of  at 
least  0.7. 

C.     THE  SYNTHESIS  MODEL 

The  synthesis  model  is  the  primary  interface  between  the  control  design 
and  the  LQR  algorithm.  At  the  heart  of  the  synthesis  model  is  a  linear  model 
of  Bluebird  developed  in  Chapter  II. 

Bluebird  has  four  control  inputs,  namely  elevator,  rudder,  ailerons,  and 
throttle.  The  elevator  and  throttle  are  natural  choices  for  controlling  x  and 
z  position  in  steady  state.  The  remaining  two  control  inputs  could  be  used 
to  control  the  lateral  variable  (y  position).  Both  rudder  and  aileron  provide 
means  of  generating  accelerations  in  the  lateral  plane.  In  fact,  rudder  is  more 
effective  at  generating  sideslip  than  aileron.  In  the  linear  plant,  lateral  position 
is  the  double  integral  of  lateral  acceleration.  Subsequently,  the  resulting  LQR 
controller  will  attempt  to  use  rudder  to  null  out  errors  in  lateral  position,  i.e., 
to  turn  the  plane.  However,  the  desired  controller  response  is  to  bank  to  turn 
using  ailerons  and  to  use  rudder  for  turn  coordination.  Furthermore,  in  the 
presence  of  wind,  it  is  desired  that  Bluebird  fly  wings  level,  crabbed  into  the 
wind,  rather  than  use  a  wing  down,  top  rudder  technique.  For  these  reasons, 
the  rudder  was  removed  as  a  control  input  to  the  linear  model. 

As  can  be  seen  from  Table  3.1,  the  dutch  roll  mode  of  Bluebird  is  lightly 
damped.  This  light  damping  of  0.111  could  pose  a  performance  problem.  Since 
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rudder  is  available  and  not  used  in  the  design  of  the  trajectory  controller,  the 
nonlinear  model  of  Bluebird  was  modified  to  include  a  yaw  damper  for  improved 
dutch  roll  damping.  Yaw  rate  was  fed  back  to  the  rudder  through  a  constant 
gain  block  with  a  value  of  0.55.  Additionally,  the  rudder  was  removed  as  an 
external  input.  Note  that  Bluebird  is  still  fully  controllable  with  the  remaining 
three  control  inputs. 

TABLE  3.1:  EIGENVALUES  OF  BLUEBIRD 


Mode 

Frequency 

Damping 

Longitudinal 

rad/sec 

Short  Period 

5.9 

0.735 

Phugoid 

0.497 

0.0344 

Lateral-Directional 

Dutch  Roll 

2.4 

0.111 

Spiral 

0.0384 

-1 

Roll  Response 

-4.572 

1 

The  nonlinear  model  of  Bluebird  with  three  inputs  and  integral  yaw  damper 
was  linearized,  as  per  Chapter  II,  returning  the  linear  model, 


s-{i  = 


=    Ax  +  Bu 


(3.7) 


where 


x  = 


and 


uvwpqrQQV    Xpos    Ypos    Zpos    6epot    6apot    6t 


iT 


pot 


W  =    [   Sdtvaior      ^aileron      ^throttle 


■\T 


This  linear  model  was  used  in  the  LQR  design.    The  eigenvalues  of  Bluebird 
with  a  yaw  damper  are  given  in  Table  3.2  where  it  can  be  seen  that  the  dutch 
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roll  mode  has  been  damped  out.  Note  that  the  number  of  states  decreases  to 
fifteen  since  the  rudder  actuator  state  was  removed. 


TABLE  3.2:  EIGENVALUES  OF  BLUEBIRD  WITH  YAW  DAMPER 


Mode 

Frequency 

Damping 

Longitudinal 

rad/sec 

Short  Period 

5.9 

0.735 

Phugoid 

0.497 

0.0344 

Lateral-Directional 

Dutch  Roll 

2.35 

0.5 

Spiral 

0.1788 

1 

Roll  Response 

-4.5686 

1 

Consider  Figure  3.3.  Here  K  is  the  controller  to  be  designed,  G  is  the 
linear  model  of  Bluebird,  and  the  block,  5,  within  the  dotted  line  is  the  synthesis 
model. 

The  signal,  iu,  represents  the  commanded  trajectory  inputs: 


w 


=    [  XpOScmd      YpOScmd      ZpOS^d      Q^nd      Qcmd      <P cmd   j 


The  signal  xi  represents  the  linear  and  angular  position  states  in  the  linear 
model. 


Xi  = 


Xpos    Ypos    Zpos    $    0    # 


iT 


The  signal  e  represents  the  errors  between  the  commanded  and  current  trajec- 
tory. The  signal  z  is  comprised  of  the  outputs  of  the  matrices,  Q,  and  R.  Since 
zero  steady  state  error  is  desired  while  tracking  a  ramp  command  in  inertia! 
position,  two  integrators  were  placed  on  each  error  signal.    This  also  ensures 
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Figure  3.3:  Synthesis  and  Analysis  Model 

perfect  tracking  of  constant  heading  trajectories  in  the  presence  of  a  constant 
wind  disturbance.  Thus  Q  was  chosen  to  be; 

£11.       £12 

Q  = 

The  values  of  CxX  were  chosen  to  place  six  transmission  zeros  from  u  to  z  at 
appropriate  locations.  If  they  are  well  chosen,  six  poles  in  the  closed  loop  plant 
will  move  very  near  to  the  placed  transmisssion  zeros.  With  that  in  mind,  the 
transmission  zeros  were  chosen  as  appropriate  target  locations  for  the  poles 
added  by  the  addition  of  the  error  states. 

The  qxx  weightings  are  used  as  a  mechanism  for  obtaining  the  desired 
command  bandwidth.  Increasing  the  value  of  qxx  increases  the  relative  propor- 
tion of  that  error  state  in  the  regulated  output  vector  z.  The  resultant  LQR 
gain  increases  the  command  bandwidth  in  that  channel  in  order  to  move  the 
controlled  state  to  its  commanded  value  more  quickly. 
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The  matrix  R  is  a  constant  diagonal  matrix  required  to  be  full  rank.  Since 
the  plant  Q  has  three  control  inputs,  R  is  of  the  form, 

/  rn     0      0 

R=\     0     r22     0 

\     0  0        T33 

The  elements  of  R  are  used  as  a  mechanism  for  selecting  the  control  bandwidth. 
An  increase  in  rxx  increases  the  relative  proportion  of  that  actuator  energy  in  the 
regulated  output  z.  The  resultant  LQR  gain  decreases  the  control  bandwidth 
of  that  control  input. 

D.     THE  DESIGN  PROCESS 

Design  requirements  given  in  the  previous  section  are  SISO  in  nature. 
They  are  expressed  as  bandwidth  limitations  of  the  individual  actuators  and 
rise  time  and  damping  characteristics  along  the  command  channels.  Note,  the 
rise  time  is  inversely  proportional  to  the  command  bandwidth.  The  following 
LQR  design  process  provided  a  means  of  obtaining  a  multivariable  solution  to 
achieve  SISO  design  specifications. 

With  an  appropriate  linear  representation  of  Bluebird  and  a  synthesis 
model  that  incorporated  well  placed  transmission  zeros,  the  design  "knobs" 
were  adjusted  in  order  to  meet  performance  requirements.  The  design  "knobs" 
are  the  elemental  weightings,  qxx  and  rxx,  in  the  Q  and  R  matrices.  The  design 
process  iterated  through  the  following  steps. 

1.  Initially  let  qxx  equal  1.  Iteratively  determine  weights  for  R  to  sat- 
isfy control  loop  bandwidth  requirements.  Increasing  rxx  decreases  the  control 
bandwidth  along  that  channel. 

2.  With  R  from  step  1,  iteratively  determine  weights  for  qxx  to  satisfy 
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command  loop  bandwidth  requirements.  Increasing  qxx  increases  the  command 
bandwidth  along  that  channel  and  decreases  the  rise  time. 

3.  If  it  is  required  to  increase  the  damping  of  a  lightly  damped  mode,  use 
an  eigenvector  decomposition  to  determine  the  primary  states  affecting  that 
mode.  Include  a  weighting  on  the  derivative  of  those  states  in  the  output  z. 

4.  Ensure  that  control-loop  bandwidths  are  still  satisfactory  with  the 
values  of  qxx  in  step  2.  It  is  possible  that  all  of  the  performance  requirements 
are  not  acheivable  within  control  bandwidth  limitations. 

5.  Connect  the  LQR  controller  to  the  linear  plant  and  evaluate  the  per- 
formance in  terms  of  command  response  and  disturbance  rejection. 

6.  Confirm  satisfaction  of  other  design  requirements,  including  damping. 

7.  If  any  step  is  unsatisfactory,  go  back  to  the  synthesis  model  and  make 
appropriate  changes.  Transmission  zeros  may  need  to  be  added,  moved,  or 
deleted.  Synthesis  model  outputs  may  need  to  be  reevaluated. 

Bode  plots  were  used  to  determine  compliance  with  the  requirements. 
After  five  iterations  through  the  seven  step  process,  the  following  values  for  Q 
and  R  matrices  resulted  in  a  controller  design  that  met  design  requirements. 


Q  = 


0.5 

0 

0  1 

0 

1 

0 

0 

0 

1 

i       02  0.01 

i       (M  0.0625 

0*4  0.0625 
a  a2 


1 


R  = 


5000  0  0 
0  1000  0 
0         0       1000 


The  transmission  zeros  created  in  the  sythesis  model  are  shown  in  Ta- 
ble 3.3. 
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TABLE  3.3:  TRANSMISSION  ZEROS  OF  SYNTHESIS  MODEL 


Channel 

Cxi 

Cx2 

Cx3 

Freqency  (rad/sec) 

Damping 

■Xpoa 

1 

0.2 

0.01 

0.1 

1 

Y 

1poa 

1 

0.4 

0.0625 

0.25 

0.8 

Zpoa 

1 

0.4 

0.0625 

0.25 

0.8 

E.     LQR  CONTROLLER  PERFORMANCE 

The  eigenvalues  of  the  feedback  interconnection  of  the  plant  and  controller 
axe  given  in  Table  3.4.  It  is  apparent  that  the  zeros  created  in  the  synthesis 
model  were  well  placed  and  attracted  the  integrators  created  by  the  addition 
of  the  error  states.  There  are  two  sets  of  lightly  damped  poles.  These  do  not 
present  a  problem  because  their  frequency  is  an  order  of  magnitude  greater 
than  the  frequency  of  the  eigenvalues  associated  with  the  trajectory  commands. 
Notice  that  the  actuator  poles  did  not  change,  indicating  that  the  control  band- 
widths  were  slower  than  the  actuator  bandwidths.  Actuator  models  provide  a 
simple  means  of  determining  if  the  control  bandwidths  exceeded  the  actuator 
bandwidths. 

Figures  3.4  through  3.6  depict  the  control-loop  bandwidths  for  the  eleva- 
tor, aileron,  and  throttle.  The  cross  coupling  between  longitudinal  and  lateral 
flight  controls  was  so  slight  that  it  is  not  shown  due  to  scale. 

Figures  3.7  through  3.9  depict  the  command-loop  bandwidth  for  step  com- 
mands in  inertia!  position.  Notice  that  there  is  some  coupling  between  X  com- 
mand and  Z  response;  the  rest  are  essentially  uncoupled. 

A  summary  of  the  resulting  command  and  control  bandwidths  achieved 
is  presented  in  Table  3.5. 
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TABLE  3.4:  EIGENVALUES  OF  THE  FEEDBACK  SYSTEM 


Mode 

Frequency  (rad/sec) 

Damping 

X  Axis  Response 

0.08 

0.92 

Y  Axis  Response 

0.21 

0.77 

Z  Axis  Response 

0.21 

0.77 

elevator 

12.1 

1.0 

aileron 

12.3 

1.0 

throttle 

12.4 

1.0 

others 

0.62 

0.79 

1.90 

0.35 

2.17 

1.0 

2.35 

0.18 

2.26 

0.68 

4.54 

1.0 

5.86 

0.74 
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Figure  3.4:  Control-Loop  Bandwidth:  Elevator  Channel 

The  response  to  two  types  of  trajectory  commands  is  of  interest.  The  first 
is  the  reponse  to  a  ramp  command  in  Y  position.  This  corresponds  to  a  change 
in  the  heading  of  the  commanded  trajectory.  The  response  in  terms  of  angle  of 
bank  and  heading  activity  is  shown  in  Figure  3.10.  The  controller  achieves  the 
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Figure  3.5:  Control-Loop  Bandwidth:  Throttle  Channel 
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Figure  3.6:  Control-Loop  Bandwidth:  Aileron  Channel 

desired  result  of  turning  Bluebird  to  the  required  heading.   The  nonminimum 
phase  response  of  the  heading  state  is  due  to  adverse  yaw. 

The  response  to  a  ramp  command  in  Z  position  corresponds  to  the  re- 
sponse to  a  change  in  flight  path  angle  of  the  commanded  trajectory.  Figure  3.11 
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Figure  3.7:  Command-Loop  Bandwidth:  X  Position  Channel 
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Figure  3.8:  Command-Loop  Bandwidth:  Y  Position  Channel 

shows  the  response  to  this  command  in  terms  of  pitch  angle,  6,  and  altitude  2. 
Note  that  the  positive  z  direction  is  down,  hence  the  negative  pitch  angle. 

Finally,  the  response  to  a  constant  wind  disturbance  is  shown  in  Fig- 
ure 3.12.  The  wind  orientation  is  such  that  it  affects  all  three  axes  of  Bluebird 
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Figure  3.9:  Command-Loop  Bandwidth:  Z  Position  Channel 


TABLE  3.5:  COMMAND  AND  CONTROL  BANDWIDTHS 


Loop 

Break  Frequency 

Control  Loop 

rad/sec 

Elevator 

4.0 

Aileron 

5.0 

Throttle 

1.0 

Command  Loop 

X  Axis 

0.75 

Y  Axis 

1.0 

Z  Axis 

1.0 

equally.  The  wind  disturbance  begins  at  t  =  0. 
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Figure  3.10:  Bank  Angle  and  Heading  Response  to  Ramp  in  Y  Command 
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Figure  3.11:  Pitch  Angle  and  Altitude  Response  to  Ramp  in  Z  Command 
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Figure  3.12:  Trajectory  Error  Due  to  a  Constant  Wind  Disturbance 
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IV.  CONTROLLER  IMPLEMENTATION 
ON  THE  NONLINEAR  PLANT 

In  Chapter  III,  a  linear  controller  was  designed  to  control  the  trajectory 
of  an  air  vehicle.  However,  this  controller  cannot  be  implemented,  as  is,  on 
the  nonlinear  plant.  It  would  only  be  effective  for  commanded  trajectories  that 
represent  relatively  small  perturbations  from  the  specific  trajectory  for  which 
it  was  designed.  Intuitively,  it  is  clear  that  the  dynamics  of  the  plant  do  not 
depend  on  the  heading  angle  ty.  Furthermore,  the  orientation  of  the  force  of 
gravity  is  the  only  change  in  the  dynamics  of  the  air  vehicle  due  to  changes  in 
$  or  0.  It  turns  out  that  these  issues  were  addressed  in  [Ref.  10],  where  a  new 
methodology  for  implementing  controllers  on  nonlinear  plants  is  proposed.  The 
method  involves  differentiating  some  of  the  inputs  to  the  controller,  hence  the 
term,  P-implementation. 

This  chapter  begins  with  a  general  description  of  the  structure  of  V- 
implementation.  Furthermore,  the  specifics  of  its  implementation  on  the  non- 
linear model  in  SIMULINK  are  discussed.  Next,  the  fidelity  of  the  nonlinear 
simulation  is  improved  by  incorporating  output  feedback.  This  step  involves 
inclusion  of  high  fidelty  sensor  models  and  Kalman  filters.  Finally,  all  of  the 
pieces  of  the  complete  nonlinear  simulation  are  brought  together  in  SIMULINK. 

A.     D-IMPLEMENTATION 

Using  the  development  in  ...Chapter  II,  the  vehicle  dynamics  can  be  ex- 
pressed in  state  space  form  as  follows: 
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jtBvcg  =  fv(Bvcg»u>B,u)+fR(A)g 
jtBu>B  =  U{Bvcg,BUB,u.) 


-  P^  ='B  R(A)Bv 


eg 


(4.1) 


where  /„  and  fw  are  continuously  differentiable  and  UG7J6  denotes  the  vector 
of  control  inputs.  To  condense  the  notation,  we  define 


Xti  .  — 


B 


'eg 


1  eg 
A 


/.(•)  := 


/.(•) 
/-(•) 


/»(•)  := 


TK)9 
0 


L(.)  := 


fl(.)      0 

0      S(.) 


where  xi  €  i?6,  x2£  B?  and  L  6  .ft6*6.  Furthermore,  we  define 


r  := 


Pc 

Ac 


€#, 


(4.2) 


Q  :=< 


as  the  vector  of  linear  and  angular  position  commands  that  Bluebird  must 

track.  With  this  notation,  the  dynamics  of  the  augmented  plant  can  be  written 

as  follows: 

*•  =  h{xv,u)  +  f2(xp) 

Xp  =  Li\Xp)Xv 

ei  =  [/  0](y2  -  r)  (.  ^ 

e2  =  [0  /](j/2  -  r)  l  '*' 

j/i  =  /i(xv,u) 

I    V2      =     Xv 

where  y\  and  y2  are  the  available  measurements,  t\  and  e2  are  the  trajectory 
errors  separated  into  linear  and  angular  components.  Notice  that  L  is  only  a 
function  of  the  orientation  vector  A  =  [0  I]xp.  For  simplicity  of  exposition,  we 
have  not  included  any  extra  dynamics  for  the  actuators  or  sensors. 
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The  set  £  of  trajectories  where  the  plant  Q  is  expected  to  operate  is  given 


by 


S    :=     {(x^Xp^ucro)  : 


^VQ     


'PO 


vo 
r0 


,r0€r}, 


/l(**0»^Po>uo)  +  /2(*po)  =  °  - 

where  T  corresponds  to  the  set  of  prescribed  linear  and  angular  positions  of 
Bluebird.  This  is  a  broad  definition  of  trimmed  flight.  Notice  that  it  does 
not  preclude  the  presence  of  an  inertia!  acceleration  due  to  centripetal  force. 
As  usual,  we  restrict  the  angular  positions  to  some  subset  of  [— 7r/2,7r/2]  x 
(— 7r/2,7r/2)  x  [— 7r/2,x/2],  as  the  inverse  of  L{.)  is  not  denned  at  0  =  ir/2. 
Notice,  from  the  definition  of  x^  G  £  and  equation  (4.1)  it  follows  that: 

Bwb  £  £  — *B  wb  =  constant 
A  G  £  — ►  A  =  constant. 

Notice  that  the  set  £  is  easily  parameterized  by  Xp„  =  tq  6  T.    Given 
(im ,  Xpo ,  u0,  r0)  €  £,  we  obtain 


yio  :=  Mzn»uo) 

^2o  :=  Spo 

el0  :=  [/  0](y2o  -  r0)  :=  0 

e2o  :=  [0  /](y2o  "  r0)  :=  0.. 

Let  6x„,  6xp,  6u,6yi,  6y2,  £r,  6ei  and  6e2  correspond  to  small  per- 
turbations of  x„,  xp,  u,  yi,  y2,  r,  ci,  and  e2  about  the  nominal  values 
x^,  Xp,,,  u0,  yi0,  y2o,  r0,  elo,  and  e2o  respectively.  The  family  of  lin- 
earized models  associated  with  the  rigid  body  Q   and  the  set  £  is  defined  as 
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Qi  :=  {Gi  (r0),  r0  <E  T},  where 

6xv  =  Ai^x,,  +  A26xp  +  B^u 

6xp  =  Z6xw  +  A4^xp 

6yi  =  C\8xv  +  i^i^u 

£/  (fo)  :=  {   £y2  =  Sxp 

fei  =  [/  0](<fy2  -  <5r) 

*e2  =  [0  TPife-Jr) 

k  <*o  =  ^(xvojrcuo), 


(4.4) 


and 


A2  = 


0   g  xf  RiAotf-^Ao) 
0  0 


A,= 


0    -'bR(Aq)Vq  x  S~l{Ao) 
0  0 


(4.5) 


Matrices  A2  and  A4  were  derived  using  the  identity  in  Appendix  C.  The  intent  of 
this  derivation  is  to  isolate  the  plant  dynamics  that  are  a  function  of  Ao.  Notice 
that  A2  represents  the  contribution  of  the  force  of  gravity  to  the  dynamics  of 
Bluebird  and  A*  represents  the  sensitivity  of  Bluebird's  trajectory  to  changes 
in  the  air  vehicle's  spatial  orientation. 

Let  r*o  €  S  be  given.  Define 


where 


A2  = 


6xv    = 

A\Sxv  +  A26xp  +  B6u 

o)=< 

6xp    = 

6yi    = 

6y2    = 
6e\     = 
6e2     = 

6xv  +  A^Sxp 

C\6xv  +  D\6u 

6xp 

[I  0](6y2-6r) 

[0  I](Sy2-6r) 

b  a0      = 

u(a\*,r0,uo), 

0    g) 
0 

<ffl(A0)" 
0 

,        A4  = 

0    -v0x 
0        0 

(4.6) 


(4.7) 


Notice  that  Equation  4.6  decribes  the  linear  model  of  Bluebird  used  for 
the  design  of  the  LQR  controller  in  the  previous  chapter,  where  ro  was  chosen 
as  the  origin  of  {/}  with  Ao  equal  to  zero.  Note,  at  this  condition,  f  R(0)  =  /. 
Recall,  the  structure  of  the  controller  developed  in  Chapter  III  is  given  by: 
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6xci    =    Bc2Sxc2 
fCi  =  \  6xc2    =    [Se,  0]T  (4.8) 

.  6u       =    Cci6xcl  +  Cc26xc2  +  Dcl6y!  +  [Da  Dc2)\6ex  6e2]T. 

Based  on  K\  ,  we  propose  to  implement  the  following  controller  for  the 
nonlinear  plant  Q  : 


fC  (A)  :=  i 


f  xel    =    Bc2L"1(A)[e1  0]T 

ic2    =    Ccix*  +  CaL-\A)[ei  0]T  +  Deiyi 
+De2L-\\)[0  e2]T 

A      =    [0/]y2 
b  u       =    Xea  +  ^oaX-^A)^!  0]T. 


(4.9) 


Figure  4.1  shows  the  block  diagram  of  fC  (A).  Notice  that  A  serves  as  a 
gain  scheduling  variable.  As  a  function  of  A,  L~l  serves  to  properly  resolve 
the  trajectory  error  at  the  input  to  the  controller.  Note,  the  controller  forms 
the  derivative  of  the  measured  outputs,  y\.  Recall,  y\  is  the  measurement  of 
the  states  xv  and  the  dynamics  of  xv  are  essentially  independent  of  the  air 
vehicle's  spatial  orientation  or  position  in  {/}.  An  integrator  at  the  output  of 
the  controller  serves  to  recover  the  properties  of  the  linear  design.  The  error  is 
formed  using  the  outputs  y2  and  the  commanded  trajectory  r.  Recall,  y2  is  the 
measurement  of  the  states  xp.  L~l  serves  to  resolve  this  error,  originally  formed 
in  {/},  in  {B}. 

It  turns  out  that  the  implementation  of  Figure  4.1  has  an  important  prop- 
erty discussed  next.  First  we  need  to  make  the  following  assumptions: 

Al.  Dim(xc2)  =  dim(u)  =  dim(y2). 

A2.  The  matrix 

si  —  Ad    Bc2 

—  Cd         CC2 

has  full  rank  at  s  =  0. 
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Figure  4.1:  Block  diagram  of  the  nonlinear  controller  K.  (A) 

A3.  The  matrix  pair  (Aci,Cci)  is  observable. 

Also  denote  by  T(Qi0  Ki  )  :  6r  — *  £y2  the  closed  loop  Unear  system  that 
results  from  connecting  Qi0  to  Ki  ,  and  by  T(Qi0  K\  )(s)  its  transfer  matrix. 
Similarly,  we  let  Ti(Q  ,  JC  )(ro)  denote  the  linearization  of  the  closed  loop  system 
T(Q  ,  K,  )  at  the  equilibrium  point  determined  by  ro  and  let  Ti(Q  ,  K  )(ro)(s) 
be  its  transfer  matrix.  Then  the  following  hold. 

•  the  feedback  systems  %{Q  ,  K  )(r0)  and  T(Qi0  ,  K,\  )  have  the  same  closed 
loop  eigenvalues; 

Ti(G  ,fC  )(r0)(s)  =       L(A0)T(glo  ,  K{  ){a)L-\ho), 

A0  =  [0  /jxpo 

Thus,  the  eigenvalues  of  the  linearizations  at  each  operating  point  axe 
preserved;  furthermore,  the  input-output  behavior  of  the  linearized  operators 
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is  preserved  in  a  well-defined  sense.  The  reader  is  referred  to  [Ref.  10]  for  a 
complete  discussion  on  approximations  to  this  method  that  avoid  using  pure 
differentiation.  The  proof  of  this  result  is  contained  in  Appendix  C. 

B.     P-EVIPLEMENTATION  OF  THE  CONTROLLER  IN  SEVIULINK 

In  general,  three  steps  are  required  to  implement  the  linear  controller 
developed  Chapter  III  on  the  nonlinear  plant.  First,  the  controller  needs  to 
form  the  derivative  of  xv  as  per  Equation  4.9.  Recall,  xv  is  equal  to  the  vec- 
tor [BVcg,BUB].  The  first  three  elements  of  the  vector,  Bvcg,  are  available  as 
processed  acceleration  outputs  from  the  IMU.  Therefore,  the  controller  need 
not  compute  the  derivative  of  BVcg  since  it  will  have  it  available  directly.  The 
derivative  of  bu>b  is  computed  by  the  controller. 

Secondly,  the  controller  needs  to  act  on  the  vectors,  ei  and  e2.  The  linear 
position  error  e\  is  formed  as  the  difference  in  commanded  and  current  posi- 
tion. The  vector  e\  is  then  multiplied  by  the  transformation  matrix  BR.  This 
effectively  resolves  the  linear  trajectory  error  in  the  body-fixed  reference  frame. 
Along  this  position  trajectory,  there  exists  a  corresponding  trajectory  of  Euler 
angles.  Recall  that  xv  =  constant  is  one  of  the  constraints  on  the  set  of  tra- 
jectories. This  includes  a  broad  range  of  flight  conditions  such  as  steady  turns, 
steady  pull-ups,  climbing  or  descending  turns,  or  constant  heading.  While  it 
is  natural  to  define  the  linear  trajectory  as  a  sequence  of  positions,  it  is  more 
convenient  to  define  the  derivatives  of  the  Euler  angles  rather  than  the  values  of 
the  angles  directly.  Consider,  for  example,  that  for  many  trajectories  of  interest 
in  £,  the  components  of  A  can  be  described  by:  6  =  0;  0  =  0;  and  $  =  desired 
turn  rate.  Furthermore,  the  relationship  between  the  rates  of  change  of  the 
Euler  angles  and  bub  is  used  as  a  means  of  resolving  the  angular  position  error 
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in  {B}.  Recall, 

bu>b  =  S-1(A)A.  (4.10) 

Therefore,  the  derivative  of  the  Euler  angle  states  is  formed  and  the  commanded 
Euler  angle  rate  is  removed  to  form  e'2.  This  error  is  resolved  in  {B}  using 
Equation  4.10.  The  integrator  at  the  end  of  the  controller  recovers  the  effect  of 
the  initial  differentiation. 

Thirdly,  the  required  error  states  are  formed  by  integrating  the  rotated 
linear  trajectory  error  vector  t\.  Figure  4.1  indicates  that  integral  action  is 
accomplished  at  the  output  of  the  controller  in  order  to  recover  the  original 
properties  of  the  linear  design.  This  accounts  for  one  of  the  integration  steps 
on  the  error.  Therefore,  only  one  additional  integrator  is  required  to  provide 
double  integration  action  on  the  trajectory  error  t\. 

Figure  4.2  shows  the  P-implemented  controller  in  SIMULINK,  file  plantl.m. 
See  Appendix  B  for  a  complete  description.  The  LQR  gain  has  been  parsed  into 
several  separate  matrices  for  clarity  of  control  action. 

C.     GENERATION  OF  THE  TRAJECTORY  COMMANDS 

The  commanded  trajectory  is  specified  with  respect  to  the  inertial  refer- 
ence frame.  At  this  point,  it  is  assumed  that  a  knowledge  of  the  air  vehicle's 
performance  capabilities  is  known  and  that  the  specified  trajectory  is  within 
those  capabilities  provided  there  is  no  wind.  The  air  vehicle  has  certain  airspeed 
restrictions  with  respect  to  the  air  mass  that  cannot  be  violated  regardless  of 
the  desired  trajectory  to  be  tracked.  These  restrictions  typically  provide  lower 
and  upper  bounds  on  the  velocity  of  the  air  vehicle  with  respect  to  the  airmass. 
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Figure  4.2:  D-Implementation  of  controller  on  Bluebird 

An  example  of  a  lower  limit  is  the  stall  speed  of  the  air  vehicle.  Such  limits  are 
usually  based  on  fundamental  physical  limitations  of  the  airframe  and  a  "fly- 
able"  trajectory  can  become  "unflyable"  under  certain  conditions.  Therefore, 
a  logic  block  positioned  between  the  commanded  trajectory  and  the  controller 
ensures  that  the  commanded  trajectory  can  be  flown  at  current  flight  conditions 
within  user  defined  indicated  airspeed  limits,  shown  in  Figure  4.3.  The  com- 
manded linear  trajectory  enters  the  block  as  a  time  stamped  position  fix  in  the 
inertial  reference  frame.  Onboard  sensors  provide  both  inertia!  velocities  from 
the  IMU  and  air  mass  velocities  from  the  pitot-static  system.  Furthermore,  a 
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dual  vane  device  instrumented  on  Bluebird  provides  both  a  and  {3  measure- 
ments. Note:  a  close  approximation  to  these  readings  could  be  obtained  from 
IMU  measurements  as  outlined  in  Chapter  II. 


Portion 

Schedule 


-J   diVdt  I 


Figure  4.3:  Commanded  Trajectory  Logic  Block 


With  these  measurements,  the  wind  vector  resolved  in  {/}  is  calculated  as: 


*W  =7n  RBVca  -'  B£,Rwv, 


eg      B  "W"     vcg 


(4.11) 


where 


w 


Vcg    = 


0 
0 


and  Vt  is  the  indicated  airspeed  obtained  from  the  pitot-static  system. 

The  commanded  indicated  airspeed  of  the  air  vehicle,  Vt-cmd,  is  calculated  as: 
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where  /ucm(f  is  the  numeric  derivative  of  the  commanded  linear  trajectory.  If  the 
commanded  indicated  airspeed  is  not  within  specified  limits,  the  commanded 
trajectory  is  altered  as  follows.  The  angles  0  and  tp  are  calculated  as  follows: 

9  =  tan-\vxlvx)  (4.12) 

and 

iJ>  =  sin-1(vy/\Ivcmd\)  (4.13) 

where  the  components  of  Ivand  are  [vx,vv,vz]T.  Note  that  the  angles  6  and  ij> 
define  the  commanded  velocity  vector's  orientation  in  {/}. 

Finally,  the  amount  that  Vt-cmd  is  outside  of  indicated  airspeed  limits  is 
subtracted  from  the  magnitude  of  Ivcmi,  resulting  in  the  magnitude  of  the  new 
commanded  inertia!  velocity,  termed  V.  V  is  then  resolved  in  {/}  according  to: 


Vmod  = 


cos(9)cos(tJ>)    —cos(0)sin(xl>)    —sin(0) 

atn(V')  cos{\j))  0 

sin(6)cos(ip)    —sin(6)sin(xj>)      cos(0) 


V  (4.14) 


where  Ivmod  is  the  new  commanded  inertia!  velocity.  This  command  is  inte- 
grated and  sent  to  the  controller  as  the  commanded  trajectory.  The  MAT  LAB 
.m  file  that  implements  this  logic  is  windlogic.m  and  can  be  found  in  Appendix 
A. 

The  net  effect  of  the  trajectory  logic  block  is  simple.  When  Bluebird  runs 
up  against  one  of  its  airspeed  limits,  the  commanded  trajectory  can  no  longer  be 
followed.  A  choice  is  made  to  maintain  the  direction  of  the  commanded  velocity 
but  change  its  magnitude.  Notice  that  this  method  does  not  affect  the  turn  rate 
associated  with  the  trajectory  and  subsequently,  no  processing  of  the  angular 
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trajectory  commands  is  required.  In  this  way,  Bluebird  is  never  commanded  to 
fly  a  trajectory  that  would  force  it  to  exceed  the  performance  limits. 

The  performance  of  the  trajectory  logic  block  can  be  seen  in  Figure  4.4. 
The  lower  limit  for  Bluebird  was  arbitrarily  choosen  as  63  feet  per  second.  Blue- 
bird is  initially  flying  due  north  at  a  ground  speed  of  73  feet  per  second,  crabbed 
into  20  feet  per  second  of  wind  from  the  west.  The  commanded  trajectory  turns 
90  degrees  to  the  east.  Notice  that  the  original  trajectory  would  result  in  an 
indicated  airspeed  of  53  fps  while  the  revised  trajectory  results  in  a  commanded 
trajectory  of  63  fps. 
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Figure  4.4:  Example  of  Commanded  Trajectory  Revision 

If  the  commanded  trajectory  is  generated  using  a  velocity  rather  than 
position  schedule,  the  differentiation  block  in  Figure  4.4  can  be  removed.  A 
velocity  schedule  is  specified  in  {/}  as  a  sine  wave  of  appropriate  magnitude, 
frequency,  and  phase  along  the  x  and  y  axis.  The  commanded  ground  speed 
corresponds  to  the  magnitude  and  the  commanded  turn  rate  corresponds  to 
the  frequency  of  the  sine  function.  Constant  heading  flight  is  a  subset  of  these 
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trajectories  where  the  frequency  is  zero.  This  method  of  generating  trajectory 
commands  makes  determining  turning  rate  (hand)  simple.  As  an  example,  con- 
sider the  case  of  generating  the  velocity  schedule  for  a  circular  flight  pattern 
at  a  ground  speed  of  100  fps  and  a  turn  rate  of  0.1  radians  per  second.  The 
derivative  of  the  commanded  trajectory,  r,  parameterized  by  time  is: 


r  :=  < 


x  axis  velocity 
y  axis  velocity 
z  axis  velocity 

6 
e 


100sin(0.1<  +  pi/2) 

100sm(0.1<  +  0) 

desired  climb  or  descent  rate 

0 

0 

0.1 


(4.15) 


The  SIMULINK  block  that  generates  the  commanded  trajectory  is  shown 
in  Figure  4.5. 
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Figure  4.5:  Generation  of  Trajectory  Commands 
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D.  STATE  FEEDBACK  TO  OUTPUT  FEEDBACK 

Up  until  this  point,  the  development  of  the  tracking  controller  has  dealt 
exclusively  with  full  state  feedback.  This  section  will  detail  a  method  whereby 
the  full  state  feedback  is  replaced  by  high  fidelity  models  of  the  onboard  sensors 
in  conjunction  with  Kalman  filters  designed  to  provide  optimal  state  estimates, 
using  onboard  sensor  data. 

1.     Sensor  Modeling 

This  work  builds  upon  sensor  models  developed  in  two  prior  theses. 
In  [Ref.  1],  a  detailed  model  of  the  inertial  measurement  unit,  IMU,  is  devel- 
oped. In  [Ref.  2],  a  detailed  model  of  the  GPS  unit  is  developed.  The  IMU  is 
a  compact,  lightweight,  low  power  unit  which  integrates  nine  sensor  measure- 
ments in  one  box.  These  sensors  are  three  axis  accelerometers,  three  axis  rate 
gyros,  pitch  and  roll  inclinometers,  and  a  magnetometer.  The  accelerometers 
are  instrument  grade,  signal  conditioned,  and  temperature  compenstated.  Full 
scale  output  is  +/  —  3  g's.  The  accelerometer's  frequency  response  is  flat  past 
100  Hz.  However,  the  antialiasing  inside  the  IMU  limits  the  effective  bandwidth 
of  all  of  the  sensors  to  20  Hz.  An  internal  initialization  program  allows  the  unit 
to  compenstate  for  accelerometer  bias  and  cross  axis  error.  Table  4.1  shows  the 
specifications  of  the  accelerometers  incorporated  into  the  sensor  model  [Ref.  8] . 

The  rate  gyros  used  by  the  IMU  are  solid  state  vibrating  element 
angular  rate  sensors.  This  relatively  new  technology  uses  no  moving  parts.  A 
piezoelectric  bender  element  is  mounted  end  to  end  but  rotated  at  a  90  degree 
angle.  The  element  fastened  to  the  base  is  resonantly  driven  such  that  the 
sense  element  swings  a  reciprocating  arc.  Under  zero  angular  rate  conditions, 
the  motion  of  the  sense  element  due  to  the  drive  element  does  not  produce  any 
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TABLE  4.1:  ACCELEROMETER  CHARACTERISTICS 


Acceleration  Range 

±2g's 

Acceleration  Bandwidth 

20  Hz 

Acceleration  Bias 

0.2%  of  Full  Scale 

Acceleration  Scale  Factor 

0.2%  of  Full  Scale 

Acceleration  Noise  Floor 

0.0005  g's 

Cross  Axis  Sensitivity 

0.5%  of  Full  Scale 

bending  of  the  sense  element.  When  a  rate  of  rotation  exists,  Coriolis  forces 
cause  momentum  to  be  transfered  into  the  plane  perpendicular  to  the  motion 
of  the  drive  element,  thus  causing  bending  of  the  sense  element.  A  pressure 
transducer  picks  up  a  signal  from  the  sense  element  when  it  is  bent  that  is 
proportional  to  the  angular  rate  with  a  phase  dependent  on  the  direction  of 
rotation.  Figure  4.6  shows  the  configuration  of  two  rate  sensors  mounted  in  a 
"tuning  fork"  configuration.  Table  4.2  shows  specifications  of  the  rate  sensors 
incorporated  into  the  sensor  model  [Ref.  8]. 


Figure  4.6:  Angular  Rate  Sensor 
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TABLE  4.2:  ANGULAR  RATE  SENSOR  CHARACTERISTICS 


Rotational  Rate  Range 

±114.6deg  /  sec 

Rotational  Rate  Bandwidth 

20  Hz 

Rotational  Rate  Bias 

2.0%  of  Full  Scale 

Rotational  Rate  Scale  Factor 

0.5%  of  Full  Scale 

Rotational  Rate  Noise  Floor 

0.05%  of  Full  Scale 

Cross  Axis  Sensitivity 

0.5%  of  Full  Scale 

The  inclinometers  utilize  a  liquid  crystal  pendulous  sensor.  It  is  a 
low  bandwidth  sensor  (  approximately  0.12  radians  per  second)  that  is  meant 
to  be  integrated  with  the  rate  sensors  for  high  bandwidth  measurements  of 
angular  position.  The  fluxgate  magnetometer  provides  heading  measurements. 
The  specifications  incorporated  into  the  Euler  angle  sensor  models  are  shown 
in  Table  4.3  [Ref.  8]. 

TABLE  4.3:    INCLINOMETER  AND  MAGNETOMETER  CHARAC- 
TERISTICS 


Pitch  and  Roll  Range 

±50  deg 

Pitch  and  Roll  Bandwidth 

1/2  Hz 

Pitch  and  Roll  Accuracy 

0.2  deg 

Heading  Range 

±180  deg 

Heading  Accuracy 

3.0  deg 

Heading  Repeatability 

0.5  deg 

Heading  Linearity 

0.5% 

GPS  provides  data  in  a  form  that  can  be  converted  to  local  tangent 
plane  position.  A  brief  summary  of  the  errors  included  in  the  GPS  model 
follows. 

GPS  Error  Sources 


Selective  Availability 


56 


-  intentional  degradation  of  pseudorange  signal  by  Departmant  of  Defense 

•  Clock  Differences 

-  drift  and  bias  in  GPS  clock 

•  Ephemeris  Error 

-  error  introduced  in  converting  pseudoranges  to  inertia!  position  fix 

Each  of  these  sensor  components  is  simulated  in  block  diagram  form  in  SIMULINK 
utilizing  internal  modeling  principles  based  on  manufacturer  specifications  and 
known  sources  of  error.  The  upper  level  SIMULINK  diagram  of  these  sensor 
models  is  shown  in  Figure  4.7  and  contained  in  the  SIMULINK  file  plant l.m. 
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Figure  4.7:  Sensor  Models  in  SIMULINK 
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2.     Kalman  Filtering 

With  the  outputs  of  the  modeled  sensors  available,  Kalman  filters 
were  developed  along  the  linear  and  angular  position  channels  to  provide  optimal 
estimates  of  the  states  for  the  controller.  An  approach  analogous  to  the  LQR 
design  was  used.  Consider  the  linear  system  described  by: 


x    =    Ax  +  Bu  +  Gw  /.  .m 

.    =    Cz  +  v  (4-16) 

where  v  and  w  are  zero  mean  white  noise  with  respective  power  spectral  densities 
of  V  and  W. 

A  gain  matrix,  L,  was  found  such  that  the  Kalman  filter  given  by: 


x    =    Ax  +  Bu  +  L(z  -  Cx)  (4.17) 

produces  an  optimal  estimate  of  x.  The  Kalman  gain  L  is  calculated  as  follows: 

L  =  YCTV~\ 
where  Y  is  positive  semidefinite  and  solves  the  algebraic  Riccati  equation: 

AY  +  YAT  -  YCTV-XCY  +  GWGT  =  0 

A  synthesis  model  was  formed  that  included  the  dynamics  of  the 
original  plant.  The  IMU  used  in  Bluebird  incorporates  an  initialization  routine 
that  removes  steady  state  bias  from  the  sensors.  Therefore,  extra  dynamics  were 
not  required  in  the  Kalman  filter  to  compensate  for  bias.  The  design  process  was 
primarily  driven  by  the  bandwidth  limitations  of  the  inclinometers  and  GPS. 
The  values  of  V  and  W  were  used  as  "knobs"  in  the  iterative  design  process. 
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For  the  GPS  sensor,  a  break  frequency  of  1/2  radian  per  second  was  desired. 
For  the  Euler  angle  sensors,  a  break  frequency  of  1/10  radian  per  second  was 
desired.  It  was  also  desired  to  wash  out  the  accelerometers  and  rate  sensor 
biases  at  low  frequencies. 

The  Kalman  filter  for  the  position  estimate  blends  processed  ac- 
celerometer  outputs  with  the  GPS  position  fixes.  Figure  4.8  shows  the  frequency 
response  along  the  x  axis  of  the  Kalman  position  filter.  The  other  two  axes  have 
similar  dynamics. 
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Figure  4.8:  Frequency  Response  of  Position  Filter 

The  Kalman  filter  for  the  Euler  angles  blends  the  outputs  of  the  angu- 
lar rate  sensors,  converted  to  Euler  angle  rates,  with  Euler  angle  measurements 
from  the  inclinometers  and  magnetometer.  Figure  4.9  shows  the  frequency  re- 
soonse  of  one  channel  of  the  Euler  angle  filter.  The  remaining  two  angle  channels 
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are  similar. 
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Figure  4.0:  Frequency  Response  of  Euler  Angle  Filter 

E.     INTEGRATION  OF  THE  FULL  NONLINEAR  SIMULATION 

The  full  nonlinear  simulation  can  now  be  pieced  together.  Recall  in  Chap- 
ter II,  the  nonlinear  rigid  body  dynamics  were  implemented  in  a  SIMULINK 
block  labeled  Equations  of  Motion.  If  the  simulation  is  expanded  to  include  the 
effects  of  a  moving  airmass,  the  dynamics  of  Bluebird  can  be  simulated  at  an 
arbitrary  night  condition.  Wind  is  usually  referenced  with  respect  to  the  iner- 
tial  reference  frame,  therefore  a  SIMULINK  block,  Wind,  is  included  in  the  full 
simulation  whose  output  is  a  vector  wv  comprised  of  the  wind  velocity  resolved 
in  {/}.  Next,  the  wind  velocity  is  resolved  in  {J5}  and  added  to  the  inertia! 
velocity  of  Bluebird,  Bvcg.  The  result  is  the  velocity  of  Bluebird  with  respect 
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to  the  airmass,  resolved  in  {B}.  This  velocity,  rather  than  Bvcg,  is  used  when 
computing  the  aerodynamic  contribution  to  the  total  external  force  and  moment 
used  in  Equation  2.25.  A  more  detailed  description  of  how  this  is  accomplished 
in  the  .m  file  state-deriv.m  is  contained  in  Appendix  A. 

All  sixteen  states  are  sent  to  a  SIMULINK  block,  Sensors,  that  models 
the  avionic  sensor  suite  onboard  Bluebird.  The  output  from  these  sensors  is 
appropriately  processed  in  a  SIMULINK  block,  Kalman  filters.  The  filtered 
output  is  directed  to  the  V-Implemented  Controller  block.  The  commands  to 
the  controller  come  from  a  trajectory  block  which  uses  the  measured  outputs 
from  the  filter  block  to  process  the  commanded  trajectory.  The  controller  gen- 
erates actuator  commands  necessary  to  maintain  the  air  vehicle  on  the  com- 
manded trajectory,  thus  completing  the  loop.  The  complete  top  level  view  of 
the  SIMULINK  nonlinear  simulation  is  shown  in  Figure  4.10  and  contained  in 
the  SIMULINK  file  plantl.m. 
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Figure  4.10:  SIMULINK  Diagram  of  the  Full  Nonlinear  Simulation 
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V.  APPLICATION  TO  THE  CONTROL  OF 

BLUEBIRD 

A.     D-IMPLEMENTED  CONTROLLER  PERFORMANCE  CHAR- 
ACTERISTICS 

The  performance  characteristics  of  the  trajectory  tracking  controller  de- 
signed in  Chapters  III  and  IV  were  evaluated  using  a  two  step  process.  First,  the 
sensor  and  filter  blocks  were  removed  and  the  controller  was  connected  directly 
to  the  nonlinear  equations  of  motion  block.  Utilizing  pure  state  feedback,  Blue- 
bird was  flown  along  two  fundamentally  different  types  of  trajectories.  These 
two  trajectories  served  as  general  examples  of  the  set  of  all  trajectories  defined 
in  Chapter  IV,  Equation  4.4.  Next,  the  sensor  and  filter  blocks  were  added. 
A  general  deparure  and  arrival  trajectory,  which  is  a  combination  of  the  two 
types  of  trajectories  tested  in  the  first  step,  was  commanded  and  flown  with 
the  controller  utilizing  output  feedback.  Data  from  this  simulation  were  used 
as  input  to  a  virtual  prototype  simulation  discussed  later. 

The  dynamic  flight  simulations  were  started  using  the  same  initial  condi- 
tion. At  this  initial  condition,  Bluebird  is  aligned  with  the  positive  x-axis  and 
trimmed  for  level  flight  at  73  fps.  The  positive  x  direction  is  considered  to  be 
heading  north.  The  mechanics  of  the  dynamic  simulation  use  a  right-hand  or- 
thogonal coordinate  system  described  in  Chapter  II.  As  such,  the  positive  y-axis 
is  pointing  east  and  the  positive  z-axis  is  pointing  down.  This  choice  is  con- 
venient from  a  computational  standpoint  since  it  coincides  with  the  body-fixed 
reference  frame  of  Bluebird  at  the  initial  condition  specified  above.  Typically, 
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however,  the  positive  z  axis  is  considered  to  be  pointing  up.  If  the  right-hand 
orthogonality  is  maintained,  the  positive  y-axis  would  then  point  towards  the 
west.  For  ease  of  visualization,  this  is  the  cooridinate  system  used  to  display 
the  results  of  the  simulations. 

The  simulations  are  intended  to  evaluate  the  capabilities  of  the  controller 
in  terms  of  the  nature  of  trajectories  that  can  be  followed  in  a  stationary  air  mass 
and  in  an  air  mass  moving  at  constant  velocity.  Two  basic  kinds  of  trimmed 
flight  serve  as  the  bases  for  the  test  trajectories.  Each  test  trajectory  is  flown  in 
no- wind  conditions,  and  then  again  with  the  wind  added  at  some  point  during 
the  flight. 

The  simplest  form  of  trimmed  flight  is  constant  velocity,  constant  heading. 
This  corresponds  to  a  trajectory  defined  by  a  ramp  command  in  inertial  posi- 
tion. This  was  the  basic  trajectory  that  the  controller  was  designed  to  track. 
Figure  5.1  shows  a  three  dimensional  plot  of  the  first  test  trajectory.  In  this 
case  the  trajectory  encompasses  30  seconds  of  flight  heading  north  at  73  fps 
followed  by  a  90  degree  turn  to  join  a  trajectory  heading  east  at  73  fps  while 
climbing  at  300  feet  per  minute.  On  the  second  flight,  a  wind  of  20  feet  per 
second  from  the  north  is  added  at  the  time  the  turn  is  commanded  (elapsed 
time  =  30  seconds). 

Figure  5.2  contains  the  first  four  graphs  of  flight  data.  The  first  graph 
shows  the  time  history  of  Bluebird's  distance  from  the  commanded  trajectory. 
The  next  three  graphs  show  the  time  history  of  the  Euler  angles.  Consider  the 
baseline  flight  (no- wind  data).  Bluebird  begins  the  turn  at  an  elapsed  time  of 
30  seconds  and  exits  the  turn  at  an  elapsed  time  of  45  seconds.  Approximately 
30  seconds  later  the  trajectory  error  is  nearly  zero.  In  the  presence  of  20  feet 
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Figure  5.1:  Test  Trajectory  #1 

per  second  of  wind  from  the  north,  the  trajectory  error  also  goes  to  zero  in 
about  the  same  period  of  time.  The  graphs  of  the  Euler  angles  indicate  that 
Bluebird  is  flying  wings  level,  crabbed  into  the  cross  wind,  which  is  the  desired 
result.  Figure  5.3  shows  the  groundspeed  and  indicated  airspeed  during  the 
flights.  Notice  that  in  both  cases,  the  commanded  groundspeed  of  73  feet  per 
second  is  eventually  maintained.  Finally,  Figure  5.4  shows  the  time  history  of 
the  control  activity  during  the  flights. 

Trimmed  flight  does  not  necessarily  have  bu>b  equal  to  zero.  For  instance, 
any  steady  turning  manuever  fits  the  definition  of  trimmed  flight  given  in  Chap- 
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Figure  5.2:  Trajectory  #1:  Position  Error  and  Euler  Angles 

ter  IV.  Figure  5.5  shows  the  three  dimensional  plot  of  the  second  test  trajectory. 
In  this  case  the  test  trajectory  is  a  helix  flown  at  73  feet  per  second.  The  turn 
rate  is  one  revolution  per  minute  and  the  helix  angle  is  4  degrees  which  corre- 
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Figure  5.3:  Trajectory  #1:  Velocity  Data 

sponds  to  a  climb  rate  for  Bluebird  of  300  feet  per  minute.  For  this  test,  no 
indicated  airspeed  limits  were  placed  on  Bluebird. 

Consider  Figure  5.6  which  shows  the  position  error  and  Euler  angle  time 
history  for  the  helix  trajectory.  Notice  that  with  no  wind  the  controller  manuev- 
ers  Bluebird  to  join  the  commanded  trajectory  with  zero  error  in  steady  state. 
The  constant  pitch  and  bank  angles  confirm  the  steady  state  performance.  On 
the  second  flight,  a  wind  of  20  feet  per  second  from  the  east  was  added  at  the 
start  of  the  helical  trajectory  (elapsed  time  equal  to  40  seconds).  Intuitively, 
it  is  clear  that  the  bank  and  pitch  angle  must  vary  as  Bluebird  traverses  the 
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Figure  5.4:  Trajectory  #1:  Control  Activity 

helix.  To  an  observer  on  Bluebird,  the  wind,  while  constant  in  {/},  represents 
a  sinusoidal  disturbance.  This  explains  the  sinusoidal  nature  of  the  position  er- 
ror around  the  helix.  Figure  5.7  shows  the  groundspeed  and  indicated  airspeed 
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Figure  5.5:  Test  Trajectory  #2 

around  the  helix  during  the  two  flights.  For  the  no- wind  flight,  the  commanded 
groundspeed  is  maintained  while  for  the  flight  in  wind,  a  one  foot  per  second 
oscillation  is  experienced.  Finally,  Figure  5.8  shows  the  time  history  of  the 
control  activity  during  the  flights. 

Four  additional  flights  were  flown  using  the  helix  trajectory  in  an  attempt 
to  ascertain  the  sensitivity  of  the  position  error  to  changes  in  commanded  turn 
rate  and  wind  velocity.  Figure  5.9  shows  the  position  error  around  the  helix 
for  three  different  turn  rates  with  a  constant  wind  of  10  feet  per  second.  The 
dashed  line  corresponds  to  a  turn  rate  of  3  degrees  per  second  or  2  minutes 
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Figure  5.6:  Trajectory  #2:  Position  Error  and  Euler  Angles 

per  revolution;  the  solid  line  corresponds  to  a  turn  rate  of  6  degrees  per  second 
or  1  minute  per  revolution;  and  the  dash  dot  line  corresponds  to  a  turn  rate 
of  9  degrees  per  second  or  40  seconds  per  revolution.  The  error  increases  with 
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Figure  5.7:  Trajectory  #2:  Velocity  Data 

increasing  turn  rate.  The  9  degree  per  second  turn  rate  corresponds  to  a  steady 
state  angle  of  bank  of  thirty  degrees,  no  wind.  It  may  not  be  desireable  to 
command  trajectories  requiring  more  than  a  certain  angle  of  bank  and  this  may 
place  an  upper  bound  on  the  error. 

Figure  5.10  shows  the  position  error  around  the  helix  for  three  different 
wind  velocities  at  a  constant  turn  rate  of  6  degrees  per  second.  The  wind  varies 
in  velocity  from  10  to  25  feet  per  second.  The  local  maxima  values  of  the 
position  error  are  proportional  to  the  magnitude  of  the  disturbance. 
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Figure  5.8:  Trajectory  #2:  Control  Activity 

B.     AN  AIRPORT  DEPARTURE  AND  ARRIVAL  FLIGHT  SIM- 
ULATION 

In  many  cases,  the  trimmed  flight  condition  of  an  air  vehicle  changes  often. 
A  good  example  of  this  would  be  a  standard  instrument  departure  or  arrival 
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Figure  5.9:  Trajectory  #2:  Position  Error  for  Varying  Turn  Rates 

to  an  airfield.  These  trajectories  are  typically  a  combination  of  constant  radius 
turns  and  wings  level  flight,  while  often  climbing  and  descending.  Trajectory 
position  errors  become  critical  when  the  air  vehicle  is  on  final  approach  with  a 
constant  heading,  constant  velocity  trajectory. 

Consider  Figure  5.11.  If  an  airfield  is  imagined  to  be  located  at  the  origin, 
then  this  trajectory  would  be  indicative  of  a  typical  departure  followed  by  a 
typical  arrival  to  that  airfield.  The  scenario  utilizes  turning  trajectories  of  three 
different  radii  connected  by  straight  line  trajectories.  The  commanded  velocity 
is  a  constant  73  feet  per  second  throughout.    Wind  is  initially  zero.    Thirty 
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Figure  5.10:  Trajectory  #2:  Position  Error  for  Varying  Wind  Velocities 

seconds  into  the  flight,  the  wind  is  added  at  10  feet  per  second  from  the  east. 
At  90  seconds  into  the  flight  the  wind  is  increased  to  20  feet  per  second  from  the 
east.  Finally,  with  Bluebird  on  final  approach  tracking  a  4  degree  glideslope, 
the  wind  is  rapidly  shifted  90  degrees  to  the  north  and  decreased  in  magnitude 
to  5  feet  per  second. 

Figure  5.12  shows  the  time  history  of  the  position  error,  wind  velocity,  and 
Euler  angles  during  the  flight.  Figure  5.13  shows  the  time  history  of  the  control 
activity.  Note,  however,  the  relative  difficulty  of  analyzing  data  of  this  nature 
as  it  is  somewhat  difficult  to  visualize.  Flight  simulation  data  was  saved  to  a  file 
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intervals):  '6' 


y  range  (ft)  2000  x  range  (ft) 

Figure  5.11:  Departure  and  Arrival  at  an  Airfield 

and  processed  for  compatiblity  as  an  input  file  for  a  3-D  visualiztion  software 
package,  Designer's  Workbench.  A  virtual  prototype  of  Monterey  Airport  and 
Bluebird  was  developed  in  [Ref.  13].  The  simulation  was  then  run  as  a  departure 
and  arrival  to  the  virtual  prototype  airfield.  In  Designer's  Workbench,  the  flight 
can  be  viewed  from  multiple  perspectives  and  virtual  prototypes  of  standard 
cockpit  displays  further  enhance  visualization.  One  possible  result  is  captured 
on  video  tape,  [Ref.  14]. 
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Figure  5.12:  Airfield  Scenario:  Position  Error  and  Euler  Angles 
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Figure  5.13:  Airfield  Scenario:  Control  Activity 
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VI.  CONCLUSIONS  AND 
RECOMMENDATIONS 

A.     CONCLUSIONS 

Based  on  the  data  presented  in  this  thesis,  the  following  conclusions  are 
drawn. 

•  SIMULINK  provides  an  effective  environment  for  the  developement  of  non- 
linear simulations  for  air  vehicles.  As  a  result  of  this  development,  a  linear 
model  of  the  plant  at  an  arbitrary  trim  condition  is  easily  obtained  for 
design  purposes. 

•  LQR  design  techniques  utilizing  a  synthesis  model  and  weighting  "knobs" 
provide  a  straight  forward  means  of  obtaining  satisfactory  controller  gains 
for  MIMO  systems  while  meeting  design  requirements. 

•  ^-Implementation  of  the  linear  trajectory  tracking  controller  allows  the 
controller  to  operate  effectively  on  the  nonlinear  plant.  In  no- wind  flight 
conditions,  trajectories  defined  by  an  arbitrary  [vq,  cjo]  are  tracked  perfectly 
in  steady  state.  For  flight  conditions  with  wind,  rejection  of  a  constant 
wind  disturbance  is  accomplished  along  the  family  of  trajectories  defined 
by  an  arbitrary  [t>o»<*>o  =  Oj.  However,  for  turning  flight,  a  constant  wind 
disturbance  in  {/}  is  seen  as  a  sinusoidal  disturbance  in  {B}  and  a  sinu- 
soidal tracking  error  results.  For  moderate  bank  angles  and  turn  rates,  the 
errors  are  usually  small. 
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•  Preprocessing  of  the  trajectory  commands  by  an  adaptive  filter  allows  for 
steady  state  control  of  the  air  vehicle's  velocity  with  respect  to  the  air  mass 
in  the  forward  path,  thus  not  affecting  stability.  With  sufficient  margins  for 
transient  deviations  in  indicated  airspeed,  this  would  allieviate  the  major 
concern  of  stalling  the  air  vehicle  when  tracking  an  inertial  trajectory  in 
wind. 

•  When  analyzing  a  nonlinear  plant  and  controller,  test  simulations  are  vital 
and  in  some  cases  the  only  means  of  performance  evaluation.  The  three 
dimensional  plots  and  time  history  graphs  are  fine  for  simple  trajectories, 
but  are  difficult  to  analyze  for  more  complex  cases.  The  capabilities  of  a 
virtual  prototyping  software  package,  like  Designer's  Workbench,  are  im- 
pressive. The  enhanced  situational  awareness  and  visualization  capabilities 
of  watching  the  designed  controller  operate  on  a  virtual  prototype  allow 
for  a  "pilot's  perspective"  feedback  not  otherwise  attainable  on  the  desk 
top. 

B.     RECOMMENDATIONS 

Based  on  the  conclusions  presented  above  and  the  experience  of  developing 
the  simulation  package  presented  in  this  thesis,  the  following  recommendations 
are  made. 

•  W^hile  the  rigid  body  equations  of  motion  are  nonlinear  with  respect  to 
the  kinematics  involved,  they  are  completely  linear  with  respect  to  the 
stability  and  control  derivatives.  The  constant  coefficient  stability  and 
control  derivatives  could  be  replaced  by  functions  when  further  flight  data 
is  available. 
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•  A  similar  design  process  used  for  the  trajectory  controller  could  be  done 
using  Tioo  design  methodology. 

•  The  trajectory  preprocessor  could  be  used  to  convert  an  inertial  fixed  tra- 
jectory into  an  air  mass  fixed  trajectory.  This  might  have  applications 
where  the  air  vehicle's  inertial  position  is  of  secondary  importance  com- 
pared to  its  performance  with  respect  to  the  air  mass. 

•  Running  simulations  real  time  in  Desiner's  Workbench  rather  than  using 
batch  post  processed  data  would  be  the  next  logical  step.  Further  work 
might  lead  to  virtual  prototype  visualizations  based  on  real-time  simula- 
tions or  downlinks  from  actual  air  vehicles. 
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APPENDIX  A:  MATLAB  FILES 

The  SIMULINK  models  of  Bluebird  (plant Lm  &  plant2.m)  use  the  fol- 
lowing MATLAB  .m  files  as  user  defined  functions. 

STATE.DERIV.M 


X  X 

%  Function  to  calculate  derivative  of  u,v,w,p,q,r  X 

%  based  on  X 

%  1:  kinematics  X 

%  2:  gravity  X 

%  3:  propulsion  X 

%  4:  aerodynamics  X 

X  X 

y,  Variables  brought  from  workspace:  X 

X  X 

y,  x  ■  [contrl  inputs,  state  variables (1  -  9),  wind  vel]   */, 

X  -   (da,de,dr,dtrt,u,v,w,p  q,r,phi,theta  psi,  wind  xyz)X 

X  x 

%  Variables  called  from  function  "blue.data"  X 

X  x 

X  rho  ■  air  density  X 

%  b  =  wing  span  X 
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X  c  =  wing  cord  '/• 

X  s  =  wing  area  X 

X  Cfo  ■  Steady  state  force  term  */• 

X  Cfu  =  Stability  derivative  for  control  inputs   X 

X  m  =  airplane  mass  X 

X  lb  ■  inertia  tensor  matrix  (body  frame)  X 

X  To  =  Thrust  scale  term  % 

X  Pe  =  Engine  postion  matrix  X 

X  X 

X  % 

X  X 


function  accel  =  state_deriv(x) 


XXXXXX  Function  call  to  get  the  aircraft  data 


[uO,wO,rho,Cfx,Cfo,Cfu,Cfxdot,s,b,c,m,Pe,To,Ib]  =  blue.data; 


XXXXXX  seperate  the  combined  vector  into  seperate  elements 

u  =  [x(D;  x(2);  x(3)] ; 

dtrt  »  x(4); 

state  =  [x(5);  x(6) ;  x(7) ;  x(8);  x(9) ;  x(10)] ; 

lambda  =  [x(ll);  x(12);  x(13)] ; 
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wind  =  Cx(14);  x(15);  x(16)] ; 


y.'/.'/.'/.X'/.  calculate  velocity  wrt  the  airmass  and  form  state  vector 
'/.'/.y.'/.'/.y.  that  will  be  used  to  calculate  the  aerodynamic  forces/moments 

ias  =  u  +  wind 

statel  =  [ias(l)-uO;  ias(2);  ias(3)-w0;  x(8)  ;  x(9) ;  x(10)] ; 

y.y.y.y.y.'/.  calculate  total  velocity,  vt 

vt  -  (ias(l)*2  +  ias(2)*2  +  ias(3)~2)~ .5; 


milt       calculate  qbar 

qbar  =   .5*rho*(vt"2) ; 

XXXm  calculate  Ml 

Ml  =  diag([l/vt,  1/vt,  1/vt,  (.5*b)/vt,  (.5*c)/vt,  (.5*b)/vt]); 

lWm     calculate  M2 

M2  =  diag([0,  0,  (.5*c)/(vt~2) ,  0,  0,  0]); 

m%X%       calculate  Sprime 

83 


Sprime  =  diag([-l,  1,  -1,  b,  c,  b]*s); 

mm  calculate  Mu 

Mu  =  inv([eye(3)*m,zeros(3) ; zeros (3) ,Ib]) ; 

Wan       calculate  Tw2b 

alpha  =  ias(3)/vt; 

beta  =  ias(2)/vt; 

Tl  =  [cos(alpha),  0,  -sin(alpha);  0,1,0;  sin(alpha),  0,  cos(alpha)]; 

T2  -  [cos (beta),  -sin(beta) ,  0;  sin(beta) ,  cos (beta),  0;  0,0,1]; 

Tw2b  ■  [T1*T2,  zeros (3);  zeros (3),  T1*T2] ; 


y.y.y.y.y.y.    calculate  chi 

Chi  =  eye (6)  -  Mu*Tv2b*qbar*Sprime*Cf xdot*M2 ; 

'/.y.'/.'/X/.   calculate  Propulsion  matrix 

Prop  =  [  eye (3); 
0,-Pe(3),Pe(2); 
Pe(3),0,-Pe(l); 
-Pe(2),Pe(l),0]; 
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y.y.'/.y.'/.'/t   calculate  gravity  vector  and  rotation  matrix  {1}  to  {B} 

Rot  =  [1,  0,  -sin(lambda(2)); 

0,  cos(lambda(l)) ,  cos(lambda(2))*sin(lambda(l)) ; 

0,  -sin(lambda(l)),  cos(lambda(2))*cos(lambda(l))] ; 

Ru2b  =  [Rot ; zeros (3) ] ; 

g  =  [0;  0;  32.174]; 
FgU  =  m*g; 

'/.'/.'/. '/.'///.   put  the  components  due  to  gravity;  thrust;  and  control  surface 
'/.y.'/.'/t'/t'/.   deflections  together  for  their  contribution  to  x-dot 

thrust  =  Prop*To*dtrt ; 

gravity  =  Ru2b*FgU; 

Ctrl  ■  qbar*(Tw2b*(Sprime*(Cfo  +  (Cfu*u)))); 

xdotu  =  (Mu*(ctrl+thrust+gravity)) ; 

'/.'/.'/.y.'/.y.    calculate  kinematic  contribution 

omegax  =  [0,-state(6) ,state(5) ;state(6) ,0,-state(4) ;-state(5) ,state(4),0] ; 
¥ilb  =  (-inv(Ib))*(omegax*Ib) ; 
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Rot  ■    [-omegax,  zeros (3);  zeros (3),  wxlb] ; 


xdotrot  =  Rot*state; 


\IJ*hl%       state  vector  feedback  component  xdot 


xdotcfx  =  qbar*(Mu*(Tw2b*(Sprime*(Cfx*(Ml*statel))))); 


XXXXXX     add  three  components  of  xdot  together  and  premult  by  inv(Chi) 


xdot*  inv( Chi)* (xdot rot+xdotcfx+xdotu) ; 


XXXXXX       return  xdot 


accel«xdot ; 


EULER-B2I.M 


X  Transformation  [p  q  r]   to     lambda-dot 

X 
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function  ldot  =  eul_b2u(x) 

,/.,/.'/.'/.7.   separate  the  composite  vector  'x'  into  [p  q  r] 
%%%%%       and  [phi  theta  psi] . 

omega  =  [x(l);  x(2);  x(3)] ; 
phi=x(4); 
theta=  x(5) ; 
psi=x(6) ; 

%%%%%   calculate  the  rotation  matrix  {1}  to  {B} 
y.'/.'/.y.'/.   based  on  euler  angles 

Rb2u  =  [1,  sin(phi)*tan(theta) ,  cos  (phi)  *tan (theta)  ; 

0,  cos(phi),  -sin(phi) ; 

0,  sin(phi)*sec(theta) ,  cos (phi) *sec (theta)] ; 


VkkhU       calculate  lamda-dot 


ldot  =  Rb2u*omega; 


EULERJ2B.M 


xmxxxmxxxxmxxmmmxmmmxxmmxmmmmmm 
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y.  % 

X  Transformation  lambda-dot     to      [p  q  r]  X 

mmmmmmmmmmmmmmmmmitxmmmm 


function  omegadot  =  eul_i2b(x) 

y.'/.'/X/.   seperate  the  composite  vector  'x'  into  lambda-dot 
%*/,%%%       and  [phi  theta  psi]  . 

ldot  =  [x(l);  x(2);  x(3)] ; 

phi=x(4); 

theta=  x(5) ; 

psi»x(6); 

%%*/*%*/*         calculate  the  rotation  matrix  {B}  to  {1} 
y.y.X'/.     based  on  euler  angles 

Rb2i  =  [1,  sin(phi)*tan(theta) ,  cos  (phi)  *tan (theta)  ; 

0,  cos(phi),  -sin(phi) ; 

0,  sin(phi)*sec(theta),  cos (phi) *sec (theta)] ; 


y.'/.y.'/.    calculate  lamda-dot 

omegadot  =  inv(Rb2i)*ldot; 
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P0SITI0N-B2I.M 


mmmmmmxmmxmxxmxmmxmmxmmmmxm 

x  x 

%  From  the  workspace:                                X 

X  X 

X  1:  free  vector  'u*  resolved  in  {B}  e(l:3)       X 

X  2:  euler  angle  vector  {phi,theta,psi}  e(4:6)    X 

%  X 

%  Returns :                                          X 

X  X 

X  1:  free  vector  'u  resolved  in  {1}             X 

X  X 


function  ans  ■  pos_b2i(e) 

XXXXX   this  will  rotate  the  trajectory  error  through  phi,  theta,  psi 
XXXXX   from  {b}  to  {i}.  (3-2-1  transformation 


phi=e(4); 
theta=e(5) ; 
psi-e(6); 


m.psi  =   [cos(psi) ,sin(psi) ,0 
-sin(psi) ,cos(psi) ,0 
0,0,1]; 

m_theta  =   Ccos(theta) ,0,-sin(theta) 

0,1,0 

sin(theta) ,0, cos (that a)] ; 

m_phi  ■   [1,0,0 

0, cos (phi) , sin (phi) 

0,-sin(phi) ,cos(phi)] ; 


rotb2i  =  inv(m_phi*m_theta*m_psi) ; 


u  =   Ce(l);   e(2);   e(3)] ; 


ans  =  rotb2i*u; 


POSITIONJ2B.M 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
x  x 
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X       From  the  workspace: 

•/. 

X  1:  free  vector  'u'  resolved  in  {1}  e(l:3) 

X  2:  euler  angle  vector  {phi,theta,psi}  e(4:6) 

*/. 

'/.       Returns : 

X 

X  1:  free  vector  'u  resolved  in  {B} 

X 

mmmmmmmmxmmmmmmmitmmmmiE 

function  ans  =  pos_i2b(e) 

XXXXX   this  will  rotate  through  phi,  theta,  psi 
XXXXX   from  {i}  to  {b}.  (3-2-1  transformation) 

phi=e(4); 
theta=e(5) ; 
psi=e(6) ; 

m_psi  =  [cos(psi) ,sin(psi) ,0 
-sin(psi) ,cos(psi) ,0 
0,0,1]; 

m_theta  =  [cos (theta) ,0,-sin(theta) 
0,1,0 

91 


sin(theta) ,0,cos(theta)] ; 

m_phi  ■   [1,0,0 

0, cos (phi) , sin (phi) 

0,-sin(phi) ,cos(phi)3 ; 


roti2b  ■   (m_phi*m_theta*m_psi) ; 


u  =   Ce(l);   e(2);  e(3)]  ; 


axis  =  roti2b*u; 


LIMITXOGIC.M 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  X 

X       funtion  to  limit  trajectory  commands,  if  required  X 

X  X 

X             from  workspace:  X 

X                           1:  commanded  inertial  velocity  X 

X                           2:  inertial  wind  X 
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X  3:  lower  IAS  limit  */. 

%  4:  upper  IAS  limit  1* 

%  % 

%  returns:  revised  commanded  velocity  '/. 

x  y. 

•/•/•/•/•/•/•/•/•/•/♦/•/•/•/•/•/•/•/•/  */*/*/*/*/  WW  •/•/•/•/•/•/•/•/•/•/•/•/ •/•/»/•/•/  •/•/•/•/  WWWWWWWV'i 


function  vcom=limit(u) 

1X1*1,1     seperate  u 

vel.i  =  [u(l);u(2);u(3)]; 

wind.i  =*  [u(4),u(5),u(6)]; 

ll=u(7) ; 

ul=u(8); 

'/.'/•'/•'/•'/•   calculate  magnitud  and  direction  of  commanded  velocity 

gs=sqrt (vel.i (1)~2  +  vel_i(2)~2); 

ang=atan2(vel_i(2) ,vel_i(l)) ; 

1*1*1*1*1*1*       calculate  commanded  IAS  (steady  state) 
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vt=  sqrt((vel_i(l)+vind_i(l))~2  +(vel_i(2)+wind_i(2))~2  +  (vel_i(3)+wind_: 


y.'/.'/.'/.'/.y.    Prepare  return  variable  (may  not  be  limited) 


vcom  =  vel_i; 


*/.'/•'/•'/•%'/•    Check  limits  and  revise  if  outside 


if  vt  >  ul; 


over  =  vt  -  ul; 


vcom(l)  -  (gs  -  over)*cos(ang) ; 

vcom(2)  s  (gs  -  over)*sin(ang) ; 
end; 


if  vt  <  11; 


under  =  11  -  vt; 


vcom(l)  =  (gs  +  under) *cos (ang) ; 
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vcom(2)  =  (gs  +  under) *s in ( ang) ; 
end; 


BLUEBIRD -DATA.M 


•/•/•/•/•/•/•/•/•/•/•/  WW  •/•/•/ •/•/•/•/  WW*/  •/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/•/  •/•/•/ 


35 

X       Aircraft  data  for  Bluebird 
% 


mmmmmmmmmmummmmmmmxmxmi 


function  [uO,wO,rho,Cfx,Cfo,Cfu,Cfxdot ,s,b,c,m,Pe,To,Ib]  =  blue_data 


/././.  7. /.   trimmed  flight  speed  and  angle  of  attack 

uO  =  73.3; 

wO  =»  0; 

M.W.W.       Density:  Sea  level-  std  day 
rho  =  .0023769; 


y.V/.V/.         derivative  matrix  due  to  state  variables 
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xxxxx 


xxxxx 


rows:  [CD  CY  CL  CI  Cm  Cn] 
col:   [u  v/U  w/U  p  q  r] 


Cf x  =  [0  0  .188  0  0  0; 

0  -.31  0  0  0  .0973; 

0  0  4.22  0  3.94  0; 

0  -.0597  0  -.363  0.1; 

0  0  -1.163  0  -11.77  0; 

0  .0487  0  -.0481  0  -.0452]; 

XXXXX    derivative  matrix  due  to  control  inputs 

XXXXX    rows:   [CD  CY  CL  CI  Cm  Cn] 
XXXXX    col:   [elev  rud  ail] 

Cfu  -  [.065  0  0; 
0  .0697  0; 
.472  .0147  0; 
0  .0028  .265; 
-1.41  0  0; 
0  -.0329  -.0347]; 

XXXXX     derivative  matrix  due  to  x-dot  (alpha_dot  k   Beta.dot) 
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Cfxdot  =    [000000; 
0  0  0  0  0  0; 
0  0  1.32  0  0  0; 
0  0  0  0  0  0; 
0  0  0  0  0  0]; 

y.'/.'/.X'/.    steady  state  force  vector 

Cfo  =  [.03;  0;  .385;  0;  0;  0] ; 

XXXXX    physical  dim. 
XXXXX    WT  =55  LBS. 


m 

■ 

1. 

7095; 

s 

= 

22 

.38; 

b 

s 

12 

.42; 

c 

ss 

1. 

802; 

XXXXX    engine  data  (4  HP  motor) 

Pe  =  [0;  0;  0] ; 
To  =  [15  ;0;0]; 

y.y.'/.y.y.         inertia  tensor  matrix 

lb  «   [  10  0  0;   0   16.12     0;   0  0  7.97]; 
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APPENDIX  B:  SIMULINK  FILES 

The  nine  state  nonlinear  model  of  bluebird  is  contained  in  the  SIMULINK 
file  E0M-9.m  and  was  trimmed  at  a  flight  condition  of 

•  flight  speed  equal  to  73  fps 

•  flight  path  angle  equal  to  zero 

using  the  TRIM  command.  The  resulting  trim  values  for  the  state  vector  and 
input  vector  are: 


73.3 

0 

-0.0023 

0 

0      and  u  = 
0 
0 
0 
0 


0 

0 

0 

0.2858 


The  LINMOD  command  was  used  to  linearize  the  sixteen  state  nonlinear 
model  of  Bluebird  (contained  in  the  SIMULINK  file  EOM-16.m  and  described 
in  Chapter  III)  about  this  trim  point.  The  resulting  linear  system  is  contained 
in  the  MATLAB  file  Linear  16. mat. 

The  rudder  was  removed  as  a  control  input  (remove  the  second  column 
of  the  B  matrix)  and  the  resulting  linear  model  was  used  as  a  basis  for  the 
synthesis  model  contained  in  the  SIMULINK  file  synthesis.m.  This  synthesis 
model  was  used  to  determine  the  LQR  gain.    The  synthesis  model,  Q  and  R 
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weighting  matrices,  and  resulting  LQR  gain  is  contained  in  the  MATLAB  file 
LQR-dat.mat. 

The  full  nonlinear  simulation  is  contained  in  the  SIMULINK  file,  plantl.m. 
The  MATLAB  .m  file  simdata  loads  the  workspace  with  the  appropriate  vari- 
ables. The  file  simdata  calls  the  .m  file  trajectory. m  in  order  to  generate  the 
trajectory  schedule.  Any  changes  to  the  commanded  trajectory  or  wind  distur- 
bance schedule  can  be  made  in  trajectory,  m. 

A  version  of  the  nonlinear  simulation  that  does  not  use  the  filter  and 
sensor  blocks  is  contained  in  the  SIMULINK  file  plant2.m.  It  runs  considerably 
quicker  than  plantl.m. 
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APPENDIX  C:  ©-IMPLEMENTATION 
PROOFS  REFERENCED 

Identity.  Lei  x  (z  B?  =  const  be  given.   Then 

d 


and 


ti(IBR(A)x)  =  -IBR(A)x  x  S-\A).  (C.l) 


-£K(fR(A)x)  =  xxfR(A)S-1(\).  (C.2) 


Proof:  To  derive  both  equations  we  will  need  Poisson's  Law: 

jt(BR(A))  =B  uB  xB  R(A),  (C.3) 

and  the  following  identity: 

axb  =  -bxa  (C.4) 

for  any  vectors  a  and  6  of  compatible  dimensions.  Now,  consider 

jt(BR(A)x)  =  (jt(BR(A))x  +{,  R{A)jtx 

=B  wB  xB  R{A)x  =  -BR(A)x  xB  uB,  (C.5) 

using  equation  (C.3),  (C.4)  and  x  =  const. 
Next,  by  the  chain  rule  we  get 

i&RM*)  =  ^(^(a)x)|a 

=  ^(^(A)x)S(A)  %,.  (C.6) 

Equation  (C.l)  now  follows  by  comparing  equations  (C.5)  and  (C.6). 
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To  obtain  equation  (C.2),  consider 

j(iU(A)f R(\))  =  jt('BR(\))fR(\)  +'B  R(\)jt(fR(A))  =  0,         (C.7) 

since  BR(A)f  R(A)  =  /  VA.  Now,  using  equations  (C.3)  and  (C.7)  we  get: 

d 


dt 


(fR(A))  =  -f R(A)BuB  x 


(C.8) 


Finally,  following  the  steps  in  the  derivation  of  equation  (C.l)  we  obtain: 

^(ffi(A)x)  =  xxffl(A)5-1(A). 


Theorem  ..1  Suppose  that  assumptions  A\  through  A3  hold. 

A\.  Dim{xC2)  =  dim(u)  =  <ftm(y3). 

A2.  The  matrix 

si  -  Ad    Bc3 

—Cci         CC2 

has  full  rank  at  s  =  0. 

A3.  The  matrix  pair  (Aci,  Cc\)  is  observable. 

Then  for  each  equilibrium  point  of  Q    in  £  the  following  properties  are 
observed: 

•  the  feedback  systems  %(Q  ,  K  )(r0)  and  T(Qi0  ,  K\  )  have  the  same  closed 
loop  eigenvalues; 


Ti(Q  X  )(r0)(s)  =       L(A0)T(Glo  ,  /C,  )(s)L-\A0), 

Ao  =  [0  /]xw 
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Proof:  In  the  proof  we  set  the  controller  matrices  £)ci,  D&  D&  to  zero.  This 
does  not  change  the  results  but  considerably  simplifies  the  algebra.  Further- 
more, we  will  drop  explicit  dependence  of  the  controller  parameters  on  a.  Let 
Po 


(xvo,xPo  =      \v     ,u0,ro  €  S)  be  given.  Consider  the  feedback  interconnection 
Ao  J 

of  the  linear  plant  Qi0  (Ao)  and  linear  controller  AC/  .  The  state  matrix  F  of  this 
feedback  system  has  the  following  form: 


F:= 


A\  Ai  BjQcx  BiCci 

I  Ai       0  0 

BC\C\  0  Aci  Bc3 

C2  0         0  0 


(C.9) 


Next  we  linearize  the  feedback  interconnection  of  the  plant  Q  and  the  controller 
K  .  However,  in  order  to  that,  first  we  must  determine  the  values  of  the  con- 
troller states  xci  and  x^  along  the  trajectory  ro  £  S.  From  equation  (4.9)  we 
obtain: 


d  r,    d  r_1/k  , 

—xclo  =  AclxcU  +  Bcl—yu  +  B&L  1(A0)e0 

A 

-1( 


~nxao  =  CciaJcio  +  Cc2L    (Ao)e0 
wo  =  xc2o . 


Notice,  since  along  r0: 


eo  =  0,    t/i0  =  const,    xC20  —  u0  =  const 


we  get 


~£ci0  —  A.C\XC\0 


£ 

di 

0  =  Cc\xci0. 
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Now,  using  Assumption  A3  we  conclude  that  xc\0  =  0. 

In  order  to  compute  the  linearization  of  the  feedback  interconnection  of 
Q  and  K  ,  we  must  first  obtain  the  linearizations  of  equations  (4.3)  and  (4.9) 
about  the  operating  points  (x^^x^^uq)  G  £  and  (xci0  =  0,xc2o  =  «o,yi0  = 
0,  cq  =  y2o  —  ro  =  0)  respectively,  determined  by  ro  €  S.  The  linearization  of 
the  plant  Q  is  given  by  (4.4).  The  linearization  of  the  controller  K  has  the 
following  form: 


U    =    AelU+Bj1  +  Bc3L-1(02-p) 

£c2     =     Cclfcl  +  Cc2L~  ($2  —  p) 


(CIO) 


It  is  easy  to  verify  that  the  state  matrix  M  of  %{Q  ,  C  )(r0)  is 


M:= 


Let 


Obviously, 


M    = 


Ax 

A2 

0 

B 

L 

A|                 0 

0 

BaCiAi 

BeldA2  +  B^L-1    Ael 

BrtP\B 

0 

CcuLr            Cc\ 

0 

/    0    0    0" 

P:= 

0    L    0    0 

0    0/0 

.0    0    0/. 

• 

"/      0      0    0' 

p-'  = 

0    L-1    0    0 

0      0/0 

.000/ 

• 

)le  algebra 

it  is  easy  to  show  that 

P~lMP 

r  Ai 

A2                0 

B 

i 

A»               0 

0 

BcxCxAy 

BC\C\A2  +  Bc3    Aci    . 

BC\C\B 

0 

< 

2C2                   Cc\ 

0 

(C.ll) 
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Ai    A2 
I     Ai 

0 
0 

B 
0 

Bc\ 

d    0 
0     / 

A\    A2 
I     Aa 

+  £c3[0    /] 

Acl 

BC\C\B 

o   cc2 

ccl 

0 

The  proof  of  the  first  part  of  the  theorem  now  follows  from  Assumption  A2  and 
an  observation  that  the  matrices  F  and  M  are  in  the  form  of  the  matrices  F 
and  M.  The  proof  of  the  second  part  of  the  theorem  consists  of  the  following 
steps: 

1.  compute  the  transfer  matrix  of  the  linearization  of  the  controller  K  (A) 
using  equation  (CIO)  from  controller  inputs  6\,  02,  p  to  controller  output 

2.  apply  a  similarity  transformation 


Pi  = 


/      0 
0    L-1 


to  the  linearization  of  the  plant  Q  (=  Q\  (ro))  given  by  equation  (4.4)  and 
derive  the  transfer  matrix  from  the  control  inputs  of  the  linear  plant  rj  to 
the  outputs  6\  and  62  using  this  new  state-space  realization; 

3  compute  the  feedback  interconnection  of  the  transfer  matrices  obtained  in 
steps  1  and  2  to  get  the  final  result. 

A  simple  computation  shows  that  the  transfer  matrix  from  the  controller  inputs 
01,  #2,  p  to  the  controller  output  r/  is  given  by: 

f)(s)    =    Ccl(sl  -  Art)-\B#L-l-(b2(s)  -  p(s)) 


4-    BcA{s)  +  Cc2L-l-(b2{s)-p(s)) 

3 

Us) 
L-ip(s) 


=    Hs) 


A 


(C.12) 
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where  Ki(s)  is  the  transfer  matrix  for  the  controller  K\  . 

Applying  similarity  transformation  Pi  to  equation  (4.4)  and  computing 
the  transfer  matrix  from  77  to  0\  and  03  results  in: 


L  few 


"  Ai    A2 
I     JU 

B  • 
0 

d     0 
.  0      L 

0 

M 


(C.13) 


where  the  transfer  matrix  is  given  in  packed  matrix  form.  A  simple  observation 
shows  that 

where  G>&,(.s)  is  the  transfer  matrix  for  the  plant  Qi0  . 

Now  routine  algebra  shows  that  the  transfer  matrix  from  p  to  62  of  the 
feedback  interconnection  of  the  transfer  matrices  in  equations  (C.12)  and  (C.13) 
is  given  by: 

Ti(Q  X  )(r0)M  =  L(A0)T(Glo  ,  fC,  )(3)X"1(A0). 
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